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PREFACE 


This is the first report of the third phase of a study by 
the Hydex Corporation for the National Aeronautics and Space ° 
Administration (NASA) . The study is funded under the Conserva- 
tion and Pollution section of the AgRISTARS program, a joint 
program for Agricultural and Resources Inventory Surveys Through 
Aerospace Remote Sensing. 

The study is to develop and test technology to improve 
hydrologic modeling by incorporating remotely sensed measurements 
in operational procedures. 

The first phase of the study determined the suitability of 
present and planned remote sensing capabilities for commonly used 
hydrologic models. 

The second phase included the development of methods to use 
remotely sensed measurements with hydrologic models to ensure 
that simulated results are consistent with observed conditions. 

As part of this phase a (correlation area) method was developed 
to combine remotely sensed and standard measurements to obtain 
improved areal averaged estimates of hydrologic variables. 

The objective of this third phase is to provide user and 
programmer information to implement the techniques developed 
during Phase II. This report presents such information for the 
correlation area method (CAM) . 
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CHAPTER I 


CORRELATION AREA METHOD 


INTRODUCTION 

The correlation area method (CAM) was developed to test 
technology for improving hydrologic modeling by incorporating 
remotely sensed measurements in operational procedures. Hydex 
published four National Aeronautics and Space Administration 
(NASA) contractor reports on the results of the study: 

1. Review of Hydrologic Models for Evaluating Use of 

( 1 ) 

Remote Sensing Capabilities ; 

2. Strategies for Usincj Remotely Sensed Data in 

Hydrologic Models' ; 

3. Combining Remotely Sensed and Other Measurements 

(3) 

for Hydrologic Areal Averages ; and 

4. Creating a Bridge Between Remote Sensing and 

(4) 

Hydrologic Models 

The initial report surveyed seven commonly used hydrologic 
models and presented information on the structure, parameters, 
states, and required inputs for each model. This report 
indicates the overall complexity of each model, showing how 
runoff, soil moisture and other processes (infiltration, 
percolation, evaporation, interception, and losses) are modeled. 
This detailed knowledge of models is required to evaluate how 
additional information (in this case, remotely sensed measurement 
of hydrologic variables) might improve the models. 

Strategies to us existing and planned remotely sensed 
information in commonly used hydrologic models are in the second 
contractors report. This report analyses, the measurement 
characteristics (resolution, error, and precision in space and 
time) of six variables considered most promising for hydrologic 


1 



models. This information, coupled with the structure of the 
seven models presented in the first report, provide a sound basis 
for evaluating the usefulness of the six variables for 
operational modeling. The six variables are: 


(a) 

soil moisture; 


(b) 

impervious area; 


(c) 

land cover; 


(d) 

areal extent of 

snow cover; 

(e) 

areal extent of 

frozen ground; and 

(f) 

water equivalent 

of snow cover. 


Hydrologic variables, such as soil moisture and the amount 
of water in a snow cover (water equivalent) , are modeled over 
time by the hydrologic model. The second report suggested that 
actual measurements of hydrologic variables, such as soil 
moisture and water equivalent of the snow cover, might help keep 
the modeled states of hydrologic variables consistent with actual 
values in the real world. Remote measurements can provide 
excellent information on the areal variation of a hydrologic 
variable; whereas point insitu measurements on the ground 
generally provide a more accurate measurement at one specific 
point. Most hydrologic models model the areal average of 
hydrologic variables. Intuitively the most accurate estimate of 
the areal average can be obtained using all available 
information . 

The third report presents the fundamental concepts for the 
Correlation Area Method (CAM) for estimating areal averages of f, 

hydrologic variables from a data base with a mixture of 
measurement technologies. Since this report may not be 
available, the concepts are reviewed in the following sections of 
this chapter. 
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The results of the first phase of the study emphasized 
several reasons why remotely sensed information has not been used 
for operational hydrologic modeling. One of these is the 
d i ss imular i ty (correspondence) in the time and space averages as 
envisioned by the hydrologic model, as exist in the real world, 
and as measured by remote sensing systems. A second reason 
derives from the uncertainties in the estimates obtained from the 
actual observations and the value envisioned by the model. Thus, 
there are two sources of information, the observations and the 
model, both of which have some degree of uncertainty. 

The fourth report presented technology to bridge the gap of 
uncertainty and correspondence between remote sensing and 
hydrologic models in order to keep our models in tune with the 
real world. This is called updating a model. Figure 1-1 
demonstrates how this bridging is accomplished. The correlation 
area method is used to produce the best estimate of a hydrologic 
state (for example, the amount of soil moisture in the surface 
layer of the soil). Modifications are required in the hydrologic 
model to achieve correspondence between this real world 
measurement and the modeled states. Programs with documentation 
for updating the models of the National Weather Service River 
Forecast System (NWSRFS) will be published in a forthcoming NASA 
contract report. 

The remainder of this chapter briefly describesthe CAM 
method. For a more complete description the reader is referred 

(3) 

to the third NASA report listed above 
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BUILDING A BRIDGE 



\ 
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BRIDGING THE GAP OF UNCERTAINTY AND CORRESPONDENCE 



CRITERIA FOR METHOD 


For a 

averages of 
merits, the 


technique to be of maximum value for estimating areal 
hydrological variables from all available measure- 
following criteria should be satisfied: 


a. It should not be dependent on particular 
measurement technologies; 


b. 


It should be an objective 


technique; 


c . 


It should work regardless 
available in any one time 


of the mix of data 
step; 


d. It should explicitly recognize the sampling 
geometry of the data; 

e. It should explicitly recognize differences 
in measurement accuracy of different tech- 
nologies; and 

f. It must produce some estimate of the accuracy 
of the areal estimate. 


The correlation area method meets the criteria listed above 
for estimating areal averages from data of various sampling 
characteristics. The algorithm performs in a desirable fashion 
giving greater weight to measurements that are more accurate and 
that cover a larger portion of the basin. 


The algorithm is a heuristic procedure (an engineering 
approach) that takes liberties with a more theoretically correct 
approach for the simplicity and operational capability. 

The weight assigned to each sample in computing an areal 
average value of the hydrological variable is determined by: 

(a) the normal correlation between that type of 
measurement and the hydrological variable; and 

(b) the area of the basin that is best represented 
by that measurement. 
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Illustration of CAM Method 


To illustrate the basic principl 
There are three measurements for the 
the basin in Fig. 1-2: a point, a tr 

sample. The areal sample covers more 


es, a simple case follows, 
hydrological variable for 
ansect (line) and an areal 
than the basin of interest. 


For the purpose 
assumptions are made 
measurement with the 


of this illustration, the following 
regarding the correlation of each 
random field of a hydrological variable. 


(a) The areal measurement has a correlation of 0.7 
with the true areal average over the sample 
domain ; 

(b) The line measurement has a correlation of 0.85 
with a random point along the flight line; that 
correlation decreases on either side of the line 
as shown by the values of the parallel dashed 
lines in Fig. 1-2; and 

(c) The point measurement has a correlation of 0.94 
with the variable at the point, decreasing around 
the point as shown by the values of the circles 
in Fig. 1-2. 

The area of the basin not encompassed within the 0.7 
correlation lines for the line measurement and the 0.7 circle for 
the point measurement indicates that the areal sample provides 
the best estimate for this area. The areal sample measurements 
have a 0.7 correlation with the random field of the variable for 


this area. For the rest of the basin, the point measurement 
(within the 0.7 correlation circle) and the line measurement 
(between the 0.7. dashed lines) provide a bet.ter estimate, than _ . 

does the areal sample measurement. The areas assigned to each of 
the three measurements are shown by the shaded areas in 
Figure 1-3. 


The informational value of each measurement is dependent 
upon how much area is assigned to that measurement and on how 
well that information correlates with the random variable in that 
area. Since the correlation of the point and line measurements 
varies for different portions of these areas, the informational 
values are not constant over these areas. 
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FIG. 1-2. EXAMPLE OF CORRELATION AREA METHOD FOR 

ASSIGNING CONVENTIONAL (point) , AIRCRAFT 
(line) AND SPACE (areal) MEASUREMENTS 
TO REPRESENTATIVE BASIN AREAS. 
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The approach to determining the best average values for the 
basin includes a weighting for each area. The correlation area 
is equal to the sums of the products of the correlation 
multiplied by the portion of the area assigned to that 
measurement covered by that correlation. 

For example, the correlation area for the point sample, A , 

... . P 

in the case illustrated m Figures 1-2 and 1-3 is approximately: 

A = 0.92 X (area between the sample and 0.9 circles) 

+0.85 X (area between the 0.9 and 0.8 circles) 

+0.75 X (area between the 0.8 and 0.7 circles) 

The correlation area for the line measurement, A^ , is 
determined by the sum of the products of the average correlation 
multiplied by the area between each correlation line on either 
side of the flight line as shown in Fig. 1-2. 

In this example, the correlation area for the areal sample, 

A is 0.7 times the area of the basin assigned to the areal 
a 

measurement as shown on Fig. 1-3. 

A 

The final estimated areal average value, P, of the variable 
for the basin is simply a weighted average of the measurements. 
The weight (X.) for a measurement is equal to its correlation 
area (A.) divided by the sum of all the correlation areas for the 
basin . 



d-1) 


A measure, co , of the overall accuracy of the final estimate 
a m J 

(P) is obtained by dividing the sum of the correlation areas (A.) 

l 

by actual area of the basin (B) . 


<o 


m 


EA . 

l 

B 


( 1 - 2 ) 
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By its definition oj must be between 0 and 1 and can be 
loosely referred to as the correlation of the estimated value (P) 
with the actual average (P) of the parameter over the basin. 

The algorithm has an intuitively reasonable behavior: more 

accurate measurements get a larger sample area, a larger 
correlation area, and, thus, a larger weight than less accurate 
samples. The location of the samples also affect the shape of 
the sample area. 

When multiple line, point and areal samples are considered, 

the shape of the sample areas becomes considerably more complex 

than for the simple case shown in Fig. 1-2, but the basic concept 

is unchanged. Details on the algorithms to handle complex cases 

( 3 ) 

are described by Johnson et al 


DATA REQUIREMENTS 

Applying the correlation area method requires these 
statistical characteristics of each measurement type in order to 
properly weight it: 


(a) 

the 

correlation of 

a point measurement with 


the 

true value (Cp) 

f 


(b) 

the 

correlation of 

a transect (line) 

measure- 


ment 

with the true 

value at a point 

along the 


line 

(C %) ; 



(c) 

the 

correlation of 

an areal average 

with 


the 

true areal average (C a ); 



(d) the rate at which the point or line correlation 
decays with distance (a); and 


(e) the measure of randomness in an areal distributed 
random field ( a) . 


The sampling correlation functions listed above require that 
several covariance parameters be estimated to implement the 
correlation area technique. These covariance parameters relate 
both to the accuracy of particular measurement^ techno log i es and 
to the correlation decay in space of the random field under 
study . 
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A brief discussion pertaining to each of the data 
requirements listed above follows: 

Point , Sample Correlation, C 

P 

The value of C is the correlation of a point measurement 
p . 

with the true value at that point. Based on investigations of 

the accuracy of ground sampling technologies, it should be 

relatively easy to estimate C . Furthermore, it is expected that 

C will be generally "large", say, approximately 0.9 or better, 

P 

for most technologies. 


Line Sample Correlation , C^ 


The value of C^ is the correlation of a flight-line sample 
with the value at a point randomly located along the flight line. 
Presumably, historical ground truth experiments will be 
sufficient to estimate C^. C^ actually mixes two effects: 
first, the accuracy of the technology as it measures the true 
flight-line average; and the second the correlation of a point 
along the flight line with the flight-line average. A more 
accurate technology will increase C^, but even a perfect 
measurement technology would not guarantee a C^ of 1.0 because of 
the variability of point values along the flight line. Thus, C£ 
should be expected to be less than C . 


Areal Sample Correlation, C^ 

The value of C is the correlation of an areal sample with 
the true areal average over the sampling domain. Presumably, 
ground truth experiments will provide information to estimate C , 

<3 

but some degree of subjectivity may be needed as well. 


Decay Parameter , a 

The variable under study is assumed to be a random field. 
This random field, whether soil moisture, snow water equivalent, 
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or whatever, is assumed to have a homogeneous mean value and 
variance and isotropic covariance in space described by a simple 
exponential covariance function. The decay parameter a describes 
the value of information at one location in estimating the random 
field at another site. 

There are three possible approaches for estimating the value 
ofa: the use of historical data, real-time data, and a 

conceptual model. However, the historical approach and the 
conceptual model approach offer the most promise for operational 
estimates of a. 

Historical Data . Historical data on soil moisture or snow 
water equivalent can be analyzed to estimate a. Then the value 
of a can be assumed to apply to current conditions. The 
difficulties with this approach are: (a) procuring a historical 

data base, and (b) developing a procedure to "stratify" the data 
for different values of a related to some easy-to- identi fy 
property of current conditions. 

Real Time Data . If enough point data are available in real 
time, a value of a can be estimated for the current condition of 
the random field. Using remotely sensed data for this approach 
is rather difficult because of the sample averaging properties of 
this type of data. 

Conceptual Model . The conceptual model approach can be 
illustrated by considering soil moisture as the product of two 
random fields: the field capacity and the fraction of field 

capacity that is filled. The variability of the field capacity 
can be related to a soil map. If a conceptual model is developed 
to relate the statistics of the "fraction of field capacity" to 
some easy to estimate parameters (e.g., the antecedent 
precipitation index), it may be possible to estimate a for the 
soil moisture without actual measurements of soil moisture. 
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Measure of Randomness , o 


The same three approaches for estimating the decay 
parameter, a r can be employed to estimate the randomness measure, 
which is the standard deviation of the random field under study. 

Data Base Availability 

The data base to compute the statistical parameters for the 
correlation area method is often unavailable. However, for the 
agricultural areas of the north central plains area of the United 
States and for the central area of Russia there are sufficient 
historical data to derive initial estimates to use the 
correlation area method for soil moisture and water equivalent of 
the snow cover. Of course, the data base from remotely sensed 
measurements of hydrological variables increases there will be 
more information to estimate the required statistical parameters. 

APPLICATION OF METHOD 

The primary application for which the correlation area 
method was developed was to provide a means to integrate 
measurements from convential hydrological networks and remote 
sensing for use in operational hydrology. However, the resulting 
improved estimates of the areal averages of hydrological 
variables potentially have many other uses. For example, such 
improved estimates of areal averages of soil moisture may be used 
with large scale climatic models or agricultural models for 
predicting droughts and crop yields. 

The method can be used in a reverse sense for network design 
to determine what combined data networks would be required to 
enhance a model’s performance. For example, the accuracy of areal 
estimates of the water equivalent of the snow cover based on 
point (insitu) measurements are often poor. Having a knowledge 
of the accuracy of remote sensing techniques (i.e., aerial gamma 
radiation surveys, Carroll et al.^ or microwave, Schmugge^) 
and an analysis of the data accuracy required for a specific 
improvement in model performance, the correlation area method 



could be applied (in reverse) to define the mix of measurement 
technologies networks required to produce the desired improve- 
ment . 

DOCUMENTATION 

The remainder of this report consists of a program to 
implement the correlation area method and user and programmer 
information. 

Chapter 2 describes changes and enhancements to the CAM as 

(3) 

it was presented in the original NASA Contractor Report 

These modifications were made primarily in the interest of 

computational efficiency and to simplify the user input. Chapter 

3 describes the input and output for the program, including a 

sample run. For those interested in using the program 

information on processing input data and on the major subroutines 

see Chapter 4. Appendix A contains a complete listing of the 

FORTRAM program for implementing CAM. An errata sheet for the 

(3) 

contractor report on CAM is given in Appendix B. 
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CHAPTER II 


PROGRAM MODIFICATION 


INTRODUCTION 

In the interest of program efficiency, ease of use and 
compatibility with other programs, a number of modifications to 
the original CAM algorithm have been made to implement it as a 
computer code. This chapter describes these modifications. 

Basin Boundary Convention 

The original CAM report adopted a counterclockwise 
convention for the description of boundaries as pairs of points. 
Since the National Weather Service River Forecast System (NWSRFS) 
has adopted a clockwise convention for their basin boundaries, 
the CAM program will do likewise. This trivial change affects 
only the sign of computed areas (see Appendix I of original 
report) and makes the "forward" direction in Appendix H of the 
original report clockwise rather than counterclockwise . 


The Grid Method 


The original CAM report describes procedures to determine 
the boundary of the area within the basin associated with each 
sample which is most highly correlated with that sample--the 
"sample area". Then the appropraite sampling correlation 
function is integrated over the region within the sample area to 
produce a weighted sample area--the "correlation area"— for each 
sample. For most sample area boundaries the required integration 
cannot be performed analytically, thus the computed correlation 
area is only approximate — the degree of approximation depends on 
the accuracy of the numerical integration step. 
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An alternate method, the grid method, overlays the basin 
with a uniform square grid. Each grid point is allocated to the 
measurement with which it is most highly correlated and assigned 
a grid weight equal to the magnitude of that correlation. To 
find the basin wide areal average the sample values are weighted 
by the sum of the grid weights allocated to the sample and 
divided by the total sum of all grid weights. The accuracy 
measure is simply the sum of all grid weights divided by the 
total number of grid points on the basin. As the mesh spacing in 
the grid is decreased, the grid method will converge to the same 
results as the orginial CAM report described. 

The grid method has several advantages. It avoids, for the 
most part, the tedious computations needed to define each sample 
area boundary. It is comparatively easier to implement as a 
computer code. Its accuracy is directly related to the same 
parameter which will control its computational demands--the mesh 
size. Its numerical accuracy depends on a single parameter, the 
mesh size, rather than choices for the number of terms used to 
approximate the multitude of integrals resulting from the 
original method. (It is likely as well that one would choose for 
expediency to approximate the circular and parabolic boundaries 
of the original report with linear approximations which would 
also introduce additional numerical accuracy parameters.) 

In short, the grid method appears to be a better 
computational approach with which to implement CAM though it is 
less analytically elegant than the original sample area approach. 

Areal Sample Processing 

As envisioned in the original report, all areal-average 
samples would be processed in the same manner regardless of the 
spatial resolution of the remote sensing technology employed. 

The basin has to be divided into sub-basins each within the 
domain of a single areal sample. If two technologies were 
employed the basin would be further subdivided to account for all 
the overlapping domains of samples. These techniques work fine 
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for samples with large spatial size as envisioned in the original 
repor t--especially if only one or two technologies were employed 
at any time. But if high-resolution data and multiple 
technologies are processed in this way thousands, perhaps 
millions, of sub-basins will be created. 

To avoid the computational burden of millions of sub-basins, 
a two-step procedure will be used for processing areal sampling 
technologies. The first step is to combine all the data from 
each individual technology to estimate the areal average for the 
entire basin. The second step is to combine the estimates from 
all areal sampling technologies into a single estimate for the 
entire basin. The final result of this two-step process is 
almost identical to the original procedure and it is much more 
efficient. 

A further reduction in computational demands can be achieved 
by making a distinction between "areal samples" and "imagery". 
When the spatial resolution of the technology is small enough it 
becomes absurd to consider a pixel to be partially within the 
basin — at this point the technology is producing an imagery 
type areal samples. The "imagery" type data is input by 
specifying the center point of each pixel and a value for each 
pixel rather than locating the four corners of each sample value 
as is required for the "areal sample" type data. This 
considerably reduces input requirements for high-resolution data 
types. For the imagery type, pixels are considered to be wholly 
within or outside the basin if their center point alone is within 
or outside the basin. The two-step procedure is described in 
some detail below. 


STEP ONE: Combine All Samples From One Technology. 

The objective of this step is to combine the areal samples 

from a single technology i to find an estimate of the areal 

i 
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average over the entire basin and also estimate the correlation Ca^ 

of M a with the actual basin-wide average. There are two 
i 

alternate procedures used--one for areal samples and one for 
imagery samples as described above. 

Areal Samples ' 

First the correlation of each areal sample value must be 
adjusted for partial coverage of the basin via equation (2-1) 
below 

Cij = Cij (B n Di j )/Di j (2-1) 

where correlation of j-th areal sample of the i-th 

technology with the true areal average for the 
j-th area (probably the same for all j). 

D^j= the area (Km^) of the j-th areal sample 
of the i-th technology. 

? . . - . 

(Bf\D i j ) = the areal (Km ) within the basin and within D — 

Ci_j= adjusted correlation of the j-th sample. 

f- 


Then the areal average Ma over the basin can be estimated 

i 

from the individual sample values equation 
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Finally, the correlation of the true areal value over the 


basin with the estimate M A can be estimated via: 

i 


r 




(2-3) 


where B= basin area. 

Imagery Samples 

In fact, imagery samples can be processed in the same way 
that areal samples are processed as described above. However, a 
significant simplication results from recognizing that (for 
imagery only) : 

(800 — )= D— for j inside basin 
= 0 for j outside basin 

because the pixels of imagery data are not considered to be 
partially in the basin. This reduces equation (2-2) above to the 

form : 


Ma 

i 


2 c ij M ii 
j in J 

basin (2-4) 


T, cij 

j in 
basin 
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Further, if C^j is the same for all j pixels (which is the 
expected case) then equation (2-4) is simply 


M A 


= I E Mij 


3 

basin 


(2-5) 


where J= the number of pixels in the basin. 


Similarly equation (2-3) reduces to 


C A .= c i 

l 


( 2 - 6 ) 


where C^= C^j is the common correlation of each pixel with 
the true areal average over the domain of the 
pixel. 

STEP TWO: Combine All Technologies 

The second step in processing areal sampling technologies is 
to combine the N estimates M A . (i=l,...,N) of the basin-wide 

l 

average into a single estimate M A and also to find C , the 
associated correlation measure. By expanding the procedures de- 
scribed in the original CAM report (Appendix F) to the N 
measurement technologies the following are obtained: 


N 


m a - E b i M ; 
i=l 


(2-7) 
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1 


( 2 - 8 ) 




which completes the areal sample processing steps. 


r 
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CHAPTER III 


PROGRAM INPUT AND OUTPUT 

/ 


INTRODUCTION 

Before running the program a data file containing the basin 
data and sample values must be set up. The name of this data 
file is supplied to the program in response to a prompt. 

The input data file consists of a sequence of card images in a 
fixed field format. First, the general format for these card 
images will be described, then the field definitions for each 
data item will be given. 

General Format 


The data file card images are of two types: control cards 

and data cards. The control cards contain a charcter in column 
one, called a control code, and are otherwise blank. The control 
code tells the program what to expect on the next few data cards. 
The following is a list of control codes, with a brief 
description of their meaning: 

Control 

Code Meaning 

* The following data cards, are comments. They 

will be echo-printed, but are otherwise ignored. 
Another * control code ends the comments. 

0 The- following- data- cards- contain -comments-,-- but - 

these comments will be printed on the output 
as a title. 

B The following data cards define the basin . 

N Begin a new set of data. The following data 

cards identify the data set. 

E End of a data set. 


22 



Control 

Code 


Meaning 


S Stop. End of all data sets. 

P The following data cards contain 

L The following data cards contain 

A The following data contain areal 

I The following data cards contain 

A typical data set will look as follows: 

comment lines 

title to be printed on the output 
basin boundary data 

point data 
flight-line data 
areal data 

areal data (2nd technology) 

new title to identify second set of data 

point data 
imagery data 
areal data 

indicates end of run 

Data cards must have a blank in column 1. The data cards between 
two control cards are to be thought of as an input group . The 
format for the group is determined by the control code that 
preceeds it. The general format for each group follows. 

Basin Boundary Group : A sequency of la t i tude- long i tude 

pairs that defines the basin boundary in a clockwise 


★ 

C 

B 


P 


a) 

ui 

(0 V 
u ' 

(0 
-M 

a 


VE 
C 


<N 

a) 
w 

u < 

to 
4J 
(tf 
Q 


-N 

P 




point data. 
line data, 
da ta . 

image data. 
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manner. This input group must precede other data for 
this case. If no basin data is given, the previous basin 
data is used. 

Point Data Group ; This group of input contains all of 
the point data for a data case. Each point value is 
described by: 

- a 12-character identifier 

- latitude 

- longitude 

- sample value 

- correlation at zero distance 

- decay parameter in Km -1 

If the correlation value or decay parameter are the same 
as on the previous data card, then they need not be re- 
peated . 

Line Data Group : This input group is similar in 

format to point data. For each line sample the 
following data is supplied. 

- a 12-character identifier 

- latitude and longitude of start of line 

- latitude and longitude of end of line 

- sample value 

- zero distance correlation 

- decay parameter in Km - * 

Here too, if correlation or decay values are the same 
for several lines, these values need not be re- 
peated on the card images. 

Areal Sample Group ; Each areal sample group gives input 
for a single areal sampling technology. Sample values 
are assumed to represent nonoverlapping quardr ilateral 
areas. The technology as a whole should have a 16- 
character descriptive title and a (default) value for 
the correlation' of each areal sample with the true 
areal value. 


24 



Each individual sample then requires: 

- the sample value (numeric) 

- the correlation measure (optional) 

- (latitude, longitude) of 1st corner 

- (latitude, longitude) of 2nd corner 

- (latitude, longitude) of 3rd corner 

- (latitude, longitude) of 4th corner 

The corner points should be entered clockwise, but the 
program can reverse a counterclockwise definition. 

Imagery Sample Group: Each imagery sample group gives 

input for a single sampling technology. The first card 
image of the group should have a 16-character descrip- 
tive title and a value for the correlation of each 
pixel with the true value for the pixel. The next data 
cards define a translation table so that pixel values 
can be input as a character code instead of a numeric 
value for each pixel. The translate table consists of 
an integer, N, that gives the number of entries in the 
table, followed by character code 1, numeric value 
1, character code 2, numeric value 2, ..., character 
code N, numeric value N. The image itself is entered 
a line at a time with each line described by: 

- (latitude, longitude) of start of line 

- (latitude, longitude) of end of line 

- a character string of the indicated length giving 
values 

Specific Format (Field Definitions) 

Since data is entered in a fixed field format, it is 
important that all data appear in the correct card image columns. 
The following defines the format for each data group. The form 
of each data item is specified by showing the appropriate FORTRAN 
format description. 

Note that all latitude values are entered as decimal 

degrees, positive for. North latitude. All longitude values are 

o 

entered as decimal degrees, positive for west of 0 . The current 
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version of the program remaps all data to a polar sterographic 
grid with standard latitutde 60°N and standard latitude 105°W; 
this is appropriate for data in North America, but inadequate for 
data from other regions. The subroutine PCHRAP should be 
modified for other regions. 

Comments and Title Groups 


Anywhere in columns 2 through 80. Maximum of 50 comment 
lines . 


Basin Boundary Group 

Latitude: columns 2-9, F8.2 

Longitude: columns 10-17, F8.2 


Latitude-longitude pairs are entered one pair per card. 


Point Data Group 


Identifier: Columns 2-13, 3A4 

Latitude: Columns 14-21, F8 . 2 

Longitude: Columns 22-29, F8.2 

Sample Value: Columns 30-37, F8 . 2 
Correlation: Columns 38-44, F7.4 

Decay Para.: Columns 45-51, F7.4 

Each point is described by the information on a single 
card. If correlation or decay parameters are omitted, 
then the corresponding colums are left blank. 

Line Data Group 


Identifier: Columns 2-13, 3A4 

Lat. of Start: Columns 14-21, F8 . 2 

Long, of Start: Columns 22-29, F8.2 

Lat. of end: Columns 30-37, F8.2 

Long, of end: Columns 38-45, F8.2 

Sample value: Columns 46-53, F8 . 2 

Correlation: Columns 54-60, F7.4 

Decay Para: Columns 61-67, F7.4 
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Each line is described by the information or a single 
card. Omitted values are noted by blanks in the 
appropriate columns. 

Areal Sample Group 

Descriptive Title: Columns 2-17, 4A4 

Lat. of 1st corner: Columns 2-9, F8.2 

Long, of 1st corner: Columns 10-17, F8.2 

Lat. of 2nd corner: Columns 18-25, F8.2 

Long, of 2nd corner: Columns 26-33, F8.2 

Lat. of 3rd corner: Columns 34-41, F8.2 

Long, of 3rd corner: Columns 42-49, F8.2 

Lat. of 4th corner: Columns 50-57, F8 . 2 

Long, of 4th corner: Columns 58-65, F8.2 

Sample value: Columns 66-73, F8.2 

Correlation: Columns 74-80, F7.4 

Each quadrilateral area is defined by the information 
on a single card. If the correlation field is left 
blank, then the value from the last card on which a 
correlation value was given is used. Each group 
can have only one title, but several areal sample 
groups can be included in a data set. Thus, if areal 
data from several technologies are available, data 
from each technology can be input as a separate 
group, with an identifying label. 

Imagery Sample Group 

Description Title: Columns 2-17 

Correlation Value: Columns 18-24, F7.4 

Number of Values in Table : Columns 2-5, 14 
Character Code: Columns 6, Al 

Numeric Value: Columns 7-14, F8.2 

Charater Code: Column 15, Al 

Numeric Value: Columns 16-23, F8.2 

Character code: Column 24, Al 

Numeric Value: Columns 25-32, F8.2 

(Continue to Column 77, then begin new card with 

character code in Column 6) 
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Lat. of start: Columns 2-9, F8.2 

Long, of start: Columns 10-17, F8.2 

Lat. of end: Columns 18-25, F8 . 2 

Long, of end: Columns 26-33, F8.2 

Number of Values: Columns 34-37, 14 

Character string: Columns 38-80 

Figure 3-1 gives a sample data set. 


Program Output 


The program produces two kinds of output: a graphic 

description of the distribution of sample values in the basin and 
a table that lists, for each sampling technology, the input 
measurement value, the computed weight assigned to that value, 
the correlation coefficient, the decay parameters and a symbol 
that relates this sample to the graphic display. The basin-wide 
areal average estimate is then given, along with the overall 
accuracy estimate. This estimate will be between 0.0 and 1.0, 
with a value of 1.0 indicating an exact estimate. The sample run 
in the next section shows this output. 
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* 


THE DATA SET IS PROVIDED BY DR. E. L. PECK OF HYDEX. 
COTTONWOOD RIVER BASIN - AIRBORNE SOIL MOISTURE PROJECT 
THE BASIN HAS 9 VERTEX- POINTS. 

THERE ARE 6 POINT MEASUREMENTS, AND 16 LINE MEASUREMENTS. 


AIRBORNE SOIL MOISTURE PROJECT 

44.295 94.445 
44.168 94.964 
43.952 95.369 
44.261 95.838 
44.200 95.917 
44.235 95.995 
44.430 95.817 
44.398 95.312 
44.459 95.291 


POINT DAT AO 2 

44.296 

94.973 

27.0 

0.94 

0. 17 

POINT DATA04 

44.472 

95.488 

23.0 



POINT DATA08 

44.269 

95.967 

25.0 



POINT DATA09 

44.298 

95.648 

28.0 



POINT DATA 19 

44.111 

95.188 

23.0 



POINT DATA23 

44.308 

94.504 

26.0 



MN 501 

44.235 

95.881 

44.235 

95.959 

26.6 

MN 502 

44.239 

95.750 

44.411 

95.750 

27.1 

MN 503C 

44.244 

95.610 

44.244 

95.755 

31.5 

MN 504 

44.160 

95.619 

44.228 

95.619 

29.8 

MN 505C 

44.037 

95.451 

44.194 

95.451 

25.1 

MN 506 

44.037 

95.323 

44.037 

95.490 

21.0 

MN 507 

44.052 

95.120 

44.206 

95.120 

20.8 

MN 508C 

44.239 

95.145 

44.406 

95.145 

29.9 

MN 509 

44.393 

95.114 

44.393 

95.412 

24.6 

MN 510 

44.231 

95.271 

44.231 

95.474 

27.1 

MN 511 

44.231 

95.130 

44.231 

95.255 

28.7 

MN 512 

44.243 

94.990 

44.243 

95. 125 

22.0 

MN 513 

44.287 

94.870 

44.246 

94.990 

29.2 

MN 514C 

44.276 

94.870 

44.276 

94.760 

22.8 

MN 515 

44.206 

94.734 

44.283 

94.734 

20.0 

MN 516 

44.343 

94.516 

44.306 

94.718 

24.1 


FIGURE 3-1. A SAMPLE DATA SET 


0.17 
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A 


AREAL 

SAMPLE 







44.500 

' 96.000 

44.500 

95.250 

44.250 95.250 

44.250 

96.000 

24.000 

44.500 

95.250 

44.50 

94.500 

44.250 94.500 

44.250 

95.250 

26.000 

44.250 

96.000 

44.250 

95.250 

44.000 95.250 

44.000 

96.000 

22.000 

44.250 

95.250 

44.250 

94.500 

44.000 94.500 

44.000 

95.250 

23.000 

IMAGERY DATA 

0.25 






10A 

22.0 B 

22.50 C 

23.0 

D 23.5 E 24.0 

F 24. 

5 G 25 

.0 H : 

I 

26.0 J 

26.50 






43.750 

96.000 

44.500 

96.000 

4CDEF 




43.750 

95.500 

44.500 

95.500 

4CDEF 




43.750 

95.000 

44.500 

95.000 

4CDEF 




43.750 

94.500 

44.500 

94.500 

4CDEF 





. FIGURE 3 -1. A SAMPLE DATA SET-(Con't) 
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Sample Run 


The following sample run used the input file shown in Figure 
3-1. Note that the text enclosed between asterisks in the input 
file is printed out when the input file is read. This allows an 
immediate check to be made as to whether the correct data file 
has been read. 
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RUN PCAMGRID 


WELCOME TO CORRELATION AREA (GRID) METHOD PROGRAM 

IN ORDER TO HUN THIS PROGRAM YOU MUST ALREADY HAVE AN EXISTING 
DATA FILE CONTAINING THE BASIN DATA, POINT, LINE AND AREAL SAMPLES 

ENTER ' ' TO CONTINUE 


INPUT DATA FILE NAME: XXXXXX.DAT 
SAMPLE . DAT 

ENTER A TWO-DIGIT (01,05 - 20) OUTPUT UNIT NUMBER. 01 FOR USER TERMINAL 
01 

THE DATA SET IS PROVIDED BY DR. E. L. PECK OF HYDEX. 

COTTONWOOD RIVER BASIN - AIRBORNE SOIL MOISTURE PROJECT 
THE BASIN HAS 9 VERTEX-POINTS. 

THERE ARE 6 POINT MEASUREMENTS, AND 16 LINE MEASUREMENTS. 

WANT INPUT DATA TO BE PRINTED? ENTER Y OR N 
N 


BASIN DATA IS READ: 

PROCESSING FOR BASIN DATA COMPLETED 

BEGIN A NEW SAMPLING DATA GROUP: 

POINT SAMPLING DATA IS READ: 

PROCESSING POINT DATA COMPLETED 

LINE SAMPLING DATA IS READ: 

PROCESSING LINE DATA COMPLETED 

AREAL SAMPLING DATA IS READ: 

PROCESSING AREAL DATA COMPLETED 

IMAGERY SAMPLING DATA IS READ: 

PROCESSING IMAGERY DATA COMPLETED 

- PROCESSING FOR THIS DATA GROUP COMPLETED 

AREAL SAMPLING SUMMARY: 

SAMPLING TECH CORRELATION VALUE WEIGHT 


1 0.6417 23.81 0.8431 

2 0.2500 24.00 0.1569 

ENTER # OF GRID POINTS (INTEGER IN FOUR DIGITS) 

2500 
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SAMPLING TECH DISTRIBUTION IN THE BASIN 


*******************************$$ $9 98 ****************** ****************** 

******************************* 999999888999 $$$$$************************* 


******************************* 9999999899999 $$$$$$$$$$$$$$?************** 

*******************************999999989999$$$AAAA$$$$$$$??????????FFFF** 

***************************9999999999988899$$AAAAAA$+$$$$??????????FFFF** 

**************$$ $$$$BBB99999999999999888888$AAAAAAA+++ ????????FFFF*** 

****** 2222222$ $$$$$$$$$ 99999999999998888888 AAAAAAAA+-H- ===?????FF***** 

******222222 2$ $$$$$$$$$ 99999999999$ $8888888 AAAAAAAA+++ ====$$$******* 

***** $$22222 22 $$$$$$ $$$$$999 $$$$$$$ $8888888 AAAAAAA+++ ====$********* 

*****$$22222 2DDDDDD$$$$$$$$$$$$$$$$$8888888&AAAAAA++ =====********** 

****$$$222222DDDDDDD$$$$$$$$$$$$$$$$888888&&&&AA+++++ ====************ 

****$ $$22222 2DDDDDDD$$$$$$$$$$$$$$//////8888&&&&&&++++++ ==************** 


****$$$222222DDDDDDDD$$$$$$$000000//////#888&&&&&&&&+++$$$$$**************** 
***$$$ $22222 2DDDDDDDD$$000000000000#####8#&&&&&&&&+$ $$$****************** 
***CC$$$22222DDDDDDDD$0000000000000////M////#77&&&&&&$$$$******************* 
**CCC$$$22222DDDDD333$0000000000000////M//77777&&$$$$$********************* 
**CCC11$222223333333330000000000000#####7 77777$$$$*********************** 
*CCCC1 1 1333233333343330000000000000#?W/777 7777 $$$************************ 
*CCC1 1***3323333344444005555500000$$$EEE777777$$************************* 
CC11 1 1****222333444444$5555555$$$$$$EEEEE77 77 7 $************************** 
Cl 1 1 1*******2334444444$5555555$$$$$EEEEEEE7777*************************** 
***1 *********$$ 4444444 $ 5555555 $$ $$$EEEEEEE7 7 7 **************************** 
***************4444444$5555555$$$$$EEEEEEE77***************************** 
********* *******444444$ $5555555$$$ $EEEEEEE 7 ****************************** 
****************** 444 $$ $ 5555555 $ $$$EEEEEEE******************************* 
********************$ $$$555555566666EEEEE******************************** 
*********************$$$ 5555556666666EEE********************************* 

***********************$86555866666686$********************************** 

************************66655666666666*********************************** 

************************** 66666666666 $*********************************** 

*************************** 6666666666 ************************************ 

***************************** 666666 $************************************* 

******************************$$$$$************************************** 

********************************$$*************************************** 

# OF GRID POINTS ON THE BASIN = 1189 

THE WIDTH OF THE GRID = 0.3813 


O 
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AIRBORNE SOIL MOISTURE PROJECT 


SAMPLING 

TECH MEASUREMENT VALUE 

WEIGHT 

CORRELATION 

ALPHA 

SYMBOL 

A 


23.8399 

0.1993 

0.680 


($) 

P 


27.0000 

0.0419 

0.940 

0.170 

(A) 

P 


23.0000 

0.0023 

0.940 

0. 170 

(B) 

P 


25.0000 

0.0160 

0.940 

0.170 

(C) 

P 


28.0000 

0.0429 

0.940 

0.170 

(D) 

P 


23.0000 

0.0445 

0.940 

0.170 

(E) 

P 


26.0000 

0.0127 

0.940 

0.170 

(F) 

L 


26.6000 

0.0140 

0.850 

0.170 

(1) 

L 


27.1000 

0.0621 

0.850 

0.170 

(2) 

L 


31.5000 

0.0313 

0.850 

0.170 

(3) 

L 


29.8000 

0.0358 

0.850 

0.170 

(4) 

L 


25.1000 

0.0559 

0.850 

0.170 

(5) 

L 


21.0000 

0.0528 

0.850 

0.170 

(6) 

L 


20.8000 

0.0358 

0.850 

0.170 

(7) 

L 


29.9000 

0.0494 

0.850 

0.170 

(8) 

L 


24.6000 

0.0753 

0.850 

0.170 

(9) 

L 


27.1000 

0.0657 

0.850 

0.170 

(0) 

L 


28.7000 

0.0293 

0.850 

0.170 

(#) 

L 


22.0000 

0.0302 

0.850 

0.170 

(&) 

L 


29.2000 

0.0260 

0.850 

0.170 

(+) 

L 


22.8000 

0.0285 

0.850 

0.170 

(-) 

L 


20.0000 

0.0192 

0.850 

0.170 

(-) 

L 


24.1000 

0.0294 

0.850 

0.170 

(?) 

THE BASIN-WIDE AREAL 

AVERAGE ESTIMATE 



25.3654 



THE OVERALL 

ACCURACY 

OF THIS ESTIMATE 

= 

0.7518 




ENTER ANOTHER # OF GRID POINTS (IF NO MORE, ENTER 0) 

0 


**** STOP 


*** END *** 


OK 11:37:48 11.848 4.360 


NOTE: This run used 11.848 seconds of CPU time and 4.360 seconds of I/O time 

on a PRIME 750 computer (PRIMOS operating system). 

a 
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CHAPTER IV 


AIDS TO USERS 


PROGRAM DESCRIPTION 

This chapter describes several subtle features of the 
program. This is done solely for documentation reasons and is 
not meant to be read and understood by all users. The subtleties 
described here are important for two reasons: to simplify the 

overall program structure and to control the amount of storage 
needed for basin and sub-basin information. 

The program is organized around a set of subroutines that 
process each type of input data group. As each data group is 
read, a subroutine that can process that type of data is called. 
This subroutine, with the help of others, stores the data in one 
large array, called the boundary data array . Initially this 
array is empty. As each group of data is read, new information 
is added to the array, or previous information is combined or 
compressed. By the end of a data set, the array is again empty. 
In the next subsection, the contents of the boundary data array 
is described. Flags that are used in processing the boundary 
data and information on how data in the array is combined or 
compressed are also explained in the following sections. 


Boundary Data Array 


The boundary data array is called BDCHAN and is used to 
store information about all of the bondaries currently being 
processed. There are several different types of boundaries. 

Each boundary type has certain important information that is 
relevant to it and must be stored. The elements of BDCHAN are of 
REAL type, and are defined as follows: 

The first 7 relative positions of the boundary data array 
for each boundary hold various control information as 
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described below. The entire BDCHAN array is made up of 
multiple boundaries chained together. 


Relative 

Position Identifier Descr iption 


1 


2 


3 

4 

5 


5 


6 


TYPE 


NPTR 


MPTS 

NPTS 

AREA 


Boundary Type 

= -1.0 to indicate a boundary which 
is to be deleted. 

= 0.0 to indicate the original 
basin boundary. Only one type 
0.0 boundary will exist. 

= 1.0 to indicate a subbasin boun- 
dary. There may be multiple 
type 1.0 boundaries at any time. 

= 2.0 to indicate a special "augu- 
mented" boundary. There can be 2 
and only 2 type 2.0 boundaries 
at any time (if any) . 

Pointer to next boundary chain in 
the BDCHAN array, i.e., BDCHAN 
( I FIX (Pointer value + 0.1)) is the 
type of the next boundary data 
chain. The pointer will be set to 
-1.0 for the last boundary data 
chain. 

Maximum number of vertex points 
which this boundary can contain. 

Actual number of vertex points 
which this boundary contains. 

For Type = 0.0 or 1.0 boundaries 
only. Computed area in same units 
as x-y grid. 


NINTER For Type = 2.0 boundaries only. 

Number of intersection-type boun- 
dary points. 

CORR Computed correlation of areal 

average value to true average over 
this boundary. Set to -1.0 if not 
available. Not used for TYPE 
= 2 . 0 . 
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7 


VAL 


Computed areal average value if 
available. Not defined if CORR 

- 1 . 0 . 


The remaining relative positions hold the actual vertex 
points of the boundary. For TYPE = 0.0 or TYPE = 1.0 
each vertex point is simply an x,y pair. For a type 2.0 
boundary each vertex point has three values: x, y, and 

a "flag" which will be described in the next section. 

For type 0.0 and 1.0 boundaries, the remaining positions, 
relative to the start of this boundary data, are: 

Relative 

Position Identifier Description 


8 

X 1 

X 

position 

of 

first 

point. 

9 

Y 1 

Y 

position 

of 

first 

point . 

10 

X 2 

X 

position 

of 

second 

point 

11 

Y 2 

Y 

position 

of 

second 

point 


6+2*NPTS X NPTS x position of last point. 

7+2*NPTS Y NpTS y position of last point. 

For boundaries with Type = 2.0, the remaining positions are: 


Relative 
Pos i tion 

8 

9 

10 


Ident i f ier Description 

X ^ X position of first point. 

Y^ Y position of first point. 

Flag-^ Flag for first point (see below). 


Note that the value of the pointer, NPTR, can 
determine exactly how many locations are used 
description of each boundary. 


be used 
for the 


to 


11 

x 2 

X position 

of 

second 

poi 

12 

Y 2 

Y position 

of 

second 

poi 
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13 


Flag2 


Flag for second point. 


5 + 3 *NPTS 

6 + 3 *NPTS 
7+3*NPTS 


X NPTS 

Y NPTS 

Flag NPTS 


X position of last vertex. 

Y position of last vertex. 
Flag position of last vertex. 


Augmented Boundary Flags 


The flags are used in two different ways: one when the 

augmented boundaries are in the process of being created, and a 
second way when they are completed. When the augmented sets are 
created, each is started with a list of original boundary points, 
then augmented with all points which are the intersection points 
of two boundaries. Each intersection point appears in both 
augmented boundaries. While the augmented sets are being 
created, these intersection points are simply numbered from 1 to 
NINTER (the number of intersections) and this numeric counter is 
used as the flag value for intersections points. After the 
augmented sets are completed, these intersection point flags are 
changed into pointers to the location (1...NPTS) of the 
associated intersection point in the other augmented boundary. 


The flag value for the original points is simply a zero 
until the augmented sets are used to define subarea boundaries. 
During this process these flag values are changed from .zero (not 
used) to a -1.0 to indicate that a point has been "used" in the 
subarea definition. 


The flag values are summarized below: 

DURING AUGMENTED BOUNDARY DEFINITION 
Flag Value 
0.0 
I 

D URING SUBBASIN DEFINITION 

An original point which has been used. 

An original point which has not been used 


- 1.0 

0.0 


Description 

An original vertex point. 

The I-th intersection point, 
1 CKNINTER 
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I. ( >0 ) 


A pointer to the I-th point (1 <I<NPTS) 
in the other augmented boundary set which 
is the same intersection point as this 
one. For example. If Flag^ Q = 12 then 
the value of Flag^ in the other aug- 
mented set is 10 . 


Boundary Data Chains 

As each group of data is read, it is processed and stored in 
the array BDCHAN. The resulting "chain" of boundary data can 
have one of six forms: 

(1) Empty: 


TOP 


< empty > 


BOTTOM 


(2) Basin only: 


TOP 

(3) 


< allocated for basin boundary > 

Packed basin plus areal sample: 


BOTTOM 


TOP 

< Packed basin > 

< 4-point subbasin > 

< — empty — > 

| BOTTOM 


boundary 

for areal samples 




(4) Allocated for augmented boundaries: 


TOP 

t 

<-Packed-> 

<-4-point-> 

Allocated 

<- for augmented -> 

Allocated for 
<-augmented areal-> 

BOTTOM 

! 

basin 


basin boundary 

sample boundary 



(5) Packed Augmented Boundaries: 


TOP 

<-Packed-> 

<-4-point-> 

<-Augmented-> 

<-Augmented-> 

<- Empty- > 

BOTTOM 


Basin 


basin 

areal 




(6) Subbasin Boundaries: 


TOP j 

<-Packed-> 

<-4-point-> 

<-Augmented-> 

< -Augmented- > 

Allocated 
<- for -> 

BOTTOM 


basin 


Basin 

Areal 

Subbasins 
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The program changes the status of the boundary data chains 
in the following way: 

The program starts with BDCHAN in status (1) and returns 
to status (1) whenever a (B)asin control card is en- 
countered. After a B control card is encountered, the 
basin data is used by subroutine PCBALC to allocate the 
entire chain to a type 0 basin boundary. BDCHAN now has 
, form (2) . After the basin boundary points have all been 
processed, PCBDPK packs the basin boundary. After pack- 
ing, a 4-point type 1 boundary is set aside for any 
areal samples. This puts BDCHAN into form (3). 

As each areal sample is processed, its boundary is placed 
in the 4-point type 1 area. The empty space below this 
is allocated equally for the type 2 boundaries: The aug- 

mented basin boundary and the augmented areal sample 
boundary. This allocation requires two calls to PCBALC. 
This process gives status (4) , ready for a call to sub- 
routine PCBAUG . 

After the augmented boundaries are defined by PCBAUG, the 
array is packed into form (5) by PCBDPK. 

All of the remaining empty space is set aside for a 
single type 1 boundary to hold the subbasins as created 
by calls to PCSUBB. This is form 6. 

When all sub-basins have been found, then the augmented 
basin, augmented areal, and the allocation for sub-basin 
will all be set to type = -1 and deleted via PCBDPK. 

This restores status 3 for the next areal sample. 

Input Data Processing 

The following is a list of the subroutines that process the 
input data, with a brief description of what each does. Note 
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that the main purpose of these routines is to load data into the 
boundary data array BDCHAN. The variable NCHAN gives the number 
of elements currently in BDCHAN. 

Subroutine Name 


PCINST 

Comment Group - just echo prints until 
next * card found. 

PC INC 

Control Group. 

PCINB 

Basin Boundary Group. 

- set NCHAN= 0 

- allocate type 0 boundary 

- read boundary points 

- convert boundary from (lat, Long) 
to (x,y) 

- pack BDCHAN 

- allocate 4-point type 1 boundary 

- compute basin area (assures clock- 
wise order too) . 

PC INN 

New Data Case 

- check for valid type 0 boundary 

- read data and output options 

- initialize NPDTA, NLDTA, and 
NADTA = 0 

- set areal average and computed 
correlation of basin boundary to 
-1.0 

- set code indicating that data input 
is valid, i.e., that PCINP, PCINL, 
PCINA, and PCINI can be called. 

PCINP 

Point Data Group 

- read each value into common block 
PCPDTA and convert (lat, long) to 
(x,y). 

PCINL 

Line Data Group 

- read flight lines into common block 
PCLDTA and convert (lat, long) to 
(x,y) . 

PCINA 

Areal Data Group 

- This is a rather complex routine. 
Basically, it combines all the 
areal sample values from a single 
technology into one entry in common 
block PCADTA . 
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PCINI 


PCINE 


PC I NS 


Image Data Group 

- Combines all values from an imagery 
type technology into a single entry 
in common. 

End of Data Case 

- Perform in order: — combine all 
existing areal and image technolog- 
ies in common block PCATA into 
single basin estimate (See Sub- 
routine PCCOMB) . 

-- use grid method to combine all 
sampling technologies. See 
PCGRID. 

— output results. 

— set code indicating that data 
input is not valid until after 
"N" card. 

Stop Program (also process like end 
of data if sample values available). 
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APPENDIX A 


CAM PROGRAM LISTING 
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PROGRAM PCAMGRID 


C 

C CORRELATION AREA GRID METHOD 

C THE PURPOSE OF THIS PROGRAM IS TO OBTAIN A ESTIMATE AREAL AVERAGE 
C VALUE FOR A BOUNDED BASIN BASED ON MEASUREMENTS OBTAINED FROM 

C POINT SAMPLES, LINE SAMPLES, OR AREAL SAMPLES 

C 
C 

COMMON /PCPDTA/MPDTA , NPDTA , ALPHAP( 500 ) , CP( 500 ) , XPDTA( 500 ) , 

* YPDTA( 500 ) , CAPDTA( 500 ) , VALP( 500 ) , DESCRP( 3 , 500 ) 
COMMON /PCLDTA/MLDTA, NLDTA,ALPHAL( 100) ,DL( 100) ,XLDTA(2 , 100) , 

* YLDTA( 2 , 100) , CALDTA (100), V ALL (100),DESCRL(3,100) 
COMMON /PCADTA/MADTA , NADTA , MATDTA , NATDTA , CA( 20 ) , V ALA( 20 ) , 

* DESCRA(4 , 20) 

COMMON /PCBDCH/MCHAN,NCHAN,BDCHAN( 10000) 

C 

CHARACTER CODE* 1 , ANS 1*1, FNAME* 1 0 , TI TLE* 7 9 , INPLST* 1 
C 

IOUNIT= 1 
WRITE( 1,1) 

1 FORMAT(// IX, 'WELCOME TO CORRELATION AREA (GRID) METHOD PROGRAM’, 

* //IX, 'IN ORDER TO RUN THIS PROGRAM YOU MUST ALREADY ', 

* 'HAVE AN EXISTING' ,/ IX, 'DATA FILE CONTAINING THE BASIN ', 

* 'DATA, POINT, LINE AND AREAL SAMPLES') 

10 WRITE( 1,3) 

3 FORMAT( /IX, 'ENTER " " TO CONTINUE') 

C 

READ( 1,5) ANSI 
5 FORMAT(Al) 

IF(ANS1 .NE. ' ') GO TO 9990 
C 

40 WRITE( 1,50) 

50 FORMAT(/ IX, 'INPUT DATA FILE NAME: XXXXXX.DAT’) 

READ( 1 , 55) FNAME 
55 FORMAT(AIO) 


OPEN( UNIT=29 ,FILE=FNAME , STATUS= ’ OLD ' , ACCESS= ' SEQUENTIAL ' , 

* FORM= ' FORMATTED ' ,ERR=10,IOSTAT=IOVAL) 

C 

IF( IOVAL .EQ. 0) GO TO 70 
WRITE ( 1,60) 

60 FORMAT(/ IX, 'THERE IS A PROBLEM WITH THE DATA FILE, REENTER') 

GO TO 10 
C 

70 WRITE( 1,75) 

75 FORMAT( /IX, 'ENTER A TWO-DIGIT (01,05 - 20) OUTPUT UNIT NUMBER.', 

* ' 01 FOR USER TERMINAL ') 

READ( 1,78) IOUNIT 

78 FORMAT( 12) 

C 

100 READ( 29,5) CODE 
150 IF (CODE .NE. '*') GO TO 200 
CALL PCINST( IOUNIT, *9990) 

GO TO 100 
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c 

200 IF (CODE .NE. *C') GO TO 300 

CALL PCINC(CODE, TITLE, INPLST) 

GO TO 150 
C 

300 IF(CODE .NE. 'B') GO TO 400 
WRITE( IOUNIT,310) 

310 FORMAT(/ IX, 'BASIN DATA IS READ:') 

CALL PCINB( CODE , IOUNIT , INPLST , *9990) 

WRITE( IOUNIT ,320) 

320 FORMAT( IX, 'PROCESSING FOR BASIN DATA COMPLETED') 

GO TO 150 
C 

400 IF (CODE .NE. 'N') GO TO 500 
WRITE (IOUNIT, 4 10) 

410 FORMAT(/ IX, 'BEGIN A NEW SAMPLING DATA GROUP:’) 

CALL PCINN(CODE , IOUNIT , INPLST , *9990) 

WRITE( IOUNIT, 420) 

420 FORMAT(/ IX, 'PROCESSING FOR THIS DATA GROUP COMPLETED') 
GO TO 150 
C 

500 IF ( CODE .NE. 'E') GO TO 600 

CALL PCINE(CODE , TITLE , IOUNIT, *9990) 

GO TO 150 
C 

600 IF (CODE .NE. 'S') GO TO 1000 
CALL PCINS( IOUNIT) 

GO TO 9990 
C 

1000 WRITE( IOUNIT, 1010) 

1010 FORMAT (/ IX, 'ERROR IN CONTROL CODE') 

C 

9990 CONTINUE 
C 

STOP 

END 
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BLOCK DATA 

COMMON /PCPDTA/MPDTA , NPDTA , ALPHAP ( 500 ) , CP ( 500 ) ,XPDTA(500) , 

* YPDTA( 500 ) , CAPDTA( 500 ) , VALP( 500 ) , DESCRP( 3,500) 
COMMON /PCLDTA/MLDTA , NLDTA , ALPHAL( 100) , CL( 100) ,XLDTA( 2,100), 

* YLDTA( 2,100), CALDTA( 100),VALL(100),DESCRL(3,100) 
COMMON /PCADTA/MADTA , NADTA , MATDTA , NATDTA , CA( 20 ) , VALA( 20 ) , 

* DESCRA(4 , 20) 

COMMON /PCBDCH/MCHAN , NCHAN , BDCHAN( 10000) 


DATA MCHAN/ 10000/ 
DATA MPDTA/500/ 
DATA MLDTA/100/ 
DATA MADTA/2000/ 
DATA MATDTA/ 20/ 
DATA NGRID/400/ 

C 

END 
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C *** SUBROUNTINE PCHRAP *** 

C 

C OBJECTIVES: 

C CONVERTS THE (LATITUDE .LONGITUDE) LOCATION ON EARTH TO THE (X.Y) 
HRAP GRIP SYSTEM LOCATION. (ASSUMES NORTH POLE = (401,1601) ). 


SUBROUTINE PCHRAP (NPAIR , FLAT , FLONG , XYGRID) 

C 

DIMENSION FLAT( NPAIR) , FLONG( NPAIR) , XYGRID( 2 , NPAIR) 

C 

DATA DEGRAD/ . 0 1745329/ .EARTHR/637 1.2/, STLON/ 105 . / , STLAT/60 . / 
C 

XMESH=4.7625 
TLAT= S TLAT* DE GRAD 
RE=(EARTHR*(1.+SIN(TLAT)))/XMESH 
C 

DO 20 1=1, NPAIR 

XLAT=FLAT( I)*DEGRAD 
WLONG=(FLONG( I )+l 80 . -STLON) *DEGRAD 
R=(RE*COS(XLAT) )/ ( 1 .+SIN(XLAT) ) 

XYGRID( 1 , I)=R*SIN(WLONG)+401 . 
XYGRID(2,I)=R*COS(WLONG)+1601. 

20 CONTINUE 
C 

RETURN 

END 
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C *** SUBROUTINE P CARE A *** 

C 

C OBJECTIVES: 

C 1. FINDS THE AREA OF A GIVEN BOUNDARY 

C 2. ASSURES THAT THE BOUNDARY IS DEFINED IN A CLOCKWISE FASHION. 

C 

C PRINCIPAL VARIABLES: 


c 

NAME 

I/O 

DIMENSION 

DESCRIPTION 

c 

c 

BND 

I 

(NPV , N) 

VERTICES OF BOUNDARY 

c 


(0) 


BND( 1 , I) = X-LOCATION OF THE I-TH VERTEX 

c 




BND( 2 , I) = Y-LOCATION OF THE I-TH VERTEX 

c 




ORDER OF PTS IS REVERSED ON OUTPUT IF ISW=1 

c 

NPV 

I 

- 

NUMBEP OF DATA VALUES PER VERTEX (=2 OR 3) 

c 

N 

I 

- 

NUMBER OF VERTICES IN BOUNDARY 

c 

AREA 

0 

- 

COMPUTED AREA 

c 

ISW 

0 

- 

RETURN STATUS, =0 NORMAL RETURN 

c 




=1 PTS RETURNED IN REVERSED ORDER 

c 




=-l ERROR OCCURED 


C 

C 

SUBROUTINE PCAREA( BND , N , NPV , AREA , ISW) 

C 

DOUBLE PRECISION DAREA 
DIMENSION BND(NPV,N) 

C 

DAREA=0 . 

TEST ARGUMENTS 
ISW=-1 

IF(N .LT. 3 .OR. NPV .LT. 2) RETURN 
ISW=0 

FIND THE AREA 

DO 100 I-l.N-l 

AD=BND( 1 , 1+1 )*BND( 2 , I)-BND( 1 , I)*BND( 2,1+1) 

DAREA=DAREA+AD 
100 CONTINUE 

DARE A= DARE A+BND ( 1 , 1)*BND( 2 ,N)-BND( 1 ,N)*BND( 2,1) 

ARE A=DAREA/ 2 . 

IF (AREA .GE. 0.0) RETURN 
AREA=ABS(AREA) 

THE BOUNDARY IS IN COUNTERCLOCKWISE, ORDER OF POINTS WILL BE REVERSED 


ISW=1 

NHALF=IFIX(N/2+0. 1) 

DO 300 IUP=1 , NHALF 
IDOWN=N-IUP+l 
DO 200 1=1, NPV 
TMP=BND( I , IUP) 

BND( I , IUP)=BND( I , I DOWN) 
BND( I , IDOWN)=TMP 
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200 CONTINUE 
300 CONTINUE 

RETURN 

END 
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C *** SUBROUTINE PCIFIN *** 

C 

C OBJECTIVES: 

C DETERMINES IF A GIVEN POINT INSIDE OR OUTSIDE A BOUNDARY 
C . 

C PRINCIPAL VARIABLES: 


c 

n 

NAME 

I/O 

DIMENSION 

DESCRIPTION 

C. 

c 

X 

I 

— 

X-POSITION OF POINT TO BE TESTED 

c 

Y 

I 

- 

Y-POSITHON OF POINT TO BE TESTED 

c 

c 

c 

BND 

I 

( NPV , N) 

VERTICES OF BOUNDARY 
BND( 1,1) = X-POSITION OF THE I-TH VERTEX 
BND( 2 , I) = Y-POSITION OF THE I-TH VERTEX 

c 

NPV 

I 

- 

NUMBER OF DATA VALUES PER VERTEX (=2 OR 3) 

c 

N 

I 

- 

NUMBER OF VERTICES IN BOUNDARY 

c 

c 

INOUT 

0 


=0 FOR (X , Y) OUTSIDE OF BOUNDARY 
=1 FOR (X, Y) INSIDE BOUNDARY 


C 

C 

SUBROUTINE PCIFIN( X, Y , BND , N, NPV , INOUT) 

C 

DIMENSION BND(NPV,N) 

C 

C INSIDE OR OUTSIDE OF A BOUNDARY DEPENDS ON THE NO. OF CROSSING OF 
C ANY HALF-LINE FROM (X,Y) TO THE BOUNDARY. ODD-INSIDE, EVEN-OUTSIDE. 
C 

NCROS=0 

C 

DO 100 1=1, N 
X1=BND( 1,1) 

Y1=BND( 2,1) 

12 = 1+1 

IF( 12 .GT. N) 12=1 
X2=BND( 1,12) 

Y2=BND( 2,12) 

END POINTS OF BOUNDARY LINE FOUND: (XI ,Y1 ) ,(X2 ,Y2) 

IF(X .LT. XI .AND. X .LT. X2) GO TO 100 

IF(X .GT. XI .AND. X .GT. X2) GO TO 100 

IF( Y .LT. Y1 .AND. Y .LT. Y2) GO TO 100 

TEST X1=X2 (VERTICAL LINE) OR NOT 

IF(X1 .EQ. X2) GO TO 50 
YSTAR=(X*(Y1-Y2)+(X1*Y2-Y1*X2) )/(Xl-X2) 

IF (YSTAR .LE. Y) NCR0S=NCR0S+1 
GO TO 100 

50 IF(Y .LE. Y1 .OR. Y .LE. Y2) NCROS=NCROS+l 
100 CONTINUE 
C 

INOUT=MOD(NCROS ,2) 

C 

RETURN 

END 
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C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


*** SUBROUTINE ' PCBALC *** 

OBJECTIVES: 

ALLOCATES SPACE FOR A NEW BOUNDARY TO BDCHAN ARRAY IN COMMON BLOCK PCBDCH 


PRINCIPAL VARIABLES: 


NAME 

I/O 

DIMENSION 

DESCRIPTION 

ITYP 

I 


BOUNDARY TYPE TO BE ALLOCATED 
=0 FOR ORIGINAL BASIN BOUNDARY 
=1 FOR A SUB-BASIN BOUNDARY 
=2 FOR AN AUGMENTED BOUNDARY 

MPTS 

I/O 


ON INPUT: # OF VERTEX POINTS TO ALLOCATED 

SET TO -1 TO ALLOCATE ALL AVAILABLE SPACE 
OUTPUT: ACTUAL # OF VERTIEX POINTS ALLOCATED 

LOC 

0 

— 

LOCATION OF BOUNDARY IN BDCHAN ARRAY 
SET LOC=-l IF NO SPACE COULD BE ALLOCATED 


SUBROUTINE PCBALC( ITYP ,MPTS ,LOC) 

COMMON /PCBDCH/MCHAN , NCHAN , BDCHAN ( 10000 ) 


LOC=-l 

NOPEN=MCHAN-NCHAN 

NPV=2 

IF( ITYP .EQ. 2) NPV=3 
MAX=(NOPEN-7)/NPV 


MAX IS THE MAX # OF VERTIEX POINTS CAN BE ALLOCATED IN BDCHAN FOR THIS 
BOUNDARY. NOTE THAT ALL BOUNDARY MUST HAVE THREE OR MORE POINTS 


IF(MAX .LT. 3) RETURN 
IF(MPTS .EQ. -1) MPTS=MAX 
IF(MPTS .GT. MAX) MPTS=MAX 
LOC=NCHAN+l 

LOC IS THE FIRST AVAILABLE APACE FOR THE NEW BOUNDARY 


IF (NCHAN .EQ. 0) GO TO 100 

FIND LAST EXISTING .BOUNDARY, SET ITS POINTER TO LOC _ . 

NPTR=1 

50 LLAST=NPTR 

NPTR= IF IX ( BDCHAN ( LLAS T+ 1 ) ) 

IF(NPTR .GT. 0) GO TO 50 
BDCHAN (LLAST+1 ) =FLOAT ( LOC ) +0 . 1 

FILL VALUES FOR THIS BOUNDARY, THE SECOND POSITION IS -1 SINCE IT IS THE 
LAST 

C THE 4-TH POSITION WILL BE FILLED WHEN BOUNDARY POINTS IS READ AND SUBROUTINE 
C PCBDPK IS CALLED 
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100 BDCHAN(LOC)=FLOAT(ITYP)+O.l 
BDCHAN(LOC+l )=— 1 .0 
BDCHAN ( L0C+2 ) =FLOAT ( MPTS ) +0 . 1 
BDCHAN ( L0C+3 ) =0 . 0 
BDCHAN ( L0C+4 ) =0 . 0 
BDCHAN ( LOC-t-5 ) =- 1 . 0 
BDCHAN (L0C+6) =0.0 
C 

IF( ITYP .EQ. 2) NCHAN=NCHAN+MPTS*3+7 
IF ( ITYP .NE. 2) NCHAN=NCHAN+MPTS*2+7 
C 

RETURN 

END 




A- 9 



o o n 


C *** SUBROUTINE PCBDPK *** 

C 

C OBJECTIVES: 

C PACKS THE BOUNDARY DATA CHAINS OF ARRAY BDCHAN IN COMMON BLOCK PCBDCH 
C LEAVING ALL AVAILABLE SPACE AT THE END. BOUNDARIES WHICH HAVE TYPE=-1 

C ARE DELETED AS PART OF THE PACKING PROCESS. THE VALUE OF NCHAN IS 

C RECOMPUTED AS WELL. 

C 

C PRINCIPAL VARIABLES: 

C NAME I/O DIMENSION DESCRIPTION 

C 

C MCHAN I - DIMENSION OF ARRAY BDCHAN 

C NCHAN I/O - ACTUAL LENGTH USED IN BDCHAN 

C BDCHAN I MCHAN BOUNDARY DATA ARRAY 

C 
C 

SUBROUTINE PCBDPK 
C 

COMMON /PCBDCH/MCHAN, NCHAN, BDCHAN ( 10000) 

C 

IF( IFIX( BDCHAN( 1 ) ) .EQ. -1 .AND. BDCHAN( 2) .LT. 0.0) GO TO 140 
C 

NPTR=1 

IF( NCHAN .LE. 0) GO TO 140 
TEST NEXT BOUNDARY FOR EMPTY SPACE 
20 L0C1=NPTR 

IF (LOCI .LT. 1) GO TO 100 
NPTR=IFIX( BDCHAN(LOC 1+1 ) ) 

MPTS=IFIX(BDCHAN(LOCl+2) ) 

NPTS=IFIX( BDCHAN ( LOC 1 +3 ) ) 

ITYP=IFIX(BDCHAN(L0C1 ) ) 

IF( ITYP .LT. 0) MOVE=NPTR-LOCl 

IF( ITYP .EQ. 0 .OR. ITYP .EQ. 1) MOVE=2*(MPTS-NPTS) 

IF(ITYP .EQ. 2) MOVE=3* (MPTS-NPTS) 

IF (MOVE .LE. 0) GO TO 20 
C 

C SCROLL UP MOVE WORDS 
C FIRST MUST FIX ALL POINTERS 
C 

LOC2=LOCl 

50 IF(LOC2 .LT. 1) GO TO 55 
N=IFIX(BDCHAN(LOC2+l)) 

BDCHAN(LOC2+l )=FLOAT( IFIX(BDCHAN(LOC2+l ) )-MOVE)+0. 1 
IF(IFIX(BDCHAN(LOC2+l)) .LT. -1) BDCHAN(LOC2+l )=-l .0 
LOC2=N 
GO TO 50 
C 

55 BDCHAN(LOCl+2)=BDCHAN(LOCl+3) 

Nl=LOCl+2*MPTS+7-MOVE 

IF( ITYP .EQ. 2) N1=L0C1+3*MPTS+ 7-MOVE ~ - - - 

IF( ITYP .EQ. -1) N1=L0C1 
DO 60 I=N1 , NCHAN-MOVE 

BDCHAN( I)=BDCHAN( I+MOVE) 
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60 CONTINUE 

NCH AN=N CH AN-MOV E 
NPTR=IFIX( BDCHAN(L0C1+1 ) ) 

IF( ITYP .EQ. -I) NPTR=L0C1 
GO TO 20 
C 

100 IF(NCHAN .LE. 0) GO TO 140 
NPTR=1 

120 L0C1=NPTR 

NPTR=IFIX(BDCHAN(L0C1+1 ) ) 

IF(NPTR .GT. 0) GO TO 120 
MPTS=IFIX( BDCHAN(L0Cl+2 ) ) 

NCHAN=LOC 1+6+2 *MPTS 

IF( IFIX(BDCHAN(L0C1 )) .EQ. 2) NCHAN=NCHAN+MPTS 

IF(ITYP .LT. 0) NCHAN=L0C1-1 

RETURN 


140 NCHAN=0 
RETURN 
END 
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c *** SUBROUTINE PCINST *** 

C 

C OBJECTIVES: 

C COMMENT GROUP 

C JUST ECHO AND PRINT ALL COMMENTS UNTIL NEXT * FOUND 

C 

C 

SUBROUTINE PCINST(IOUNIT,*) 

C 

CHARACTER COMMNT*79 ,CODE*l 
C 

1 FORMAT(Al ,A79) 

2 F0RMAT(1X,A79) 

3 F0RMAT(/1X, ’TOO MANY COMMENT LINES (MAX: 50 LINES)’) 
C 

1=0 

100 1 = 1+1 

READ(29 , 1 ) CODE ,COMMNT 
IF (CODE .EQ. '*') RETURN 
WRITE( IOUNIT, 2 ) COMMNT 
IF( I .GT. 50) GO TO 9000 
C 

GO TO 100 
C 

9000 WRITE( IOUNIT, 3) 

C 

RETURN 1 
END 
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C *** SUBROUNTINE PCINB *** 

C 

C OBJECTIVES: 

C BASIN BOUNDARY GROUP 
C 1. SETS NCHAN=0 

C 2. ALLOCATES TYPE 0 BOUNDARY 

C 3. READS BASIN BOUNDARY POINTS 

C 4. CONVERTS BOUNDARY FROM (LAT..LONG.) TO (X,Y) HRAP SYSTEM 

C 5. PACKS BOUNDARY POINTS IN BDCHAN 

C 6. COMPUTES BASIN AREA 

C 7. ALLOCATES 4-POINT TYPE 1 BOUNDARY 

C 

C 

SUBROUTINE PCINB( CODE , IOUNIT , INPLST , * ) 

C 

COMMON /PCBDCH/MCHAN , NCHAN , BDCHAN( 10000) 

C 

DIMENSION FLAT( 1000 ) , FLONG( 1000) , BXY( 2,1000) 


CHARACTER CODE* 1 , INPLST* 1 
ALLOCATE TYPE 0 BOUNDARY 

NCHAN=0 /* BDCHAN IS EMPTY */ 

MPTS1=-1 

CALL PCBALC ( 0 , MPTS 1 , LOC 1 ) 

NPTS=0 

DO 100 1=1 ,MPTS1+1 

READ(29 , 1 ) CODE ,FLAT( I) , FLONG( I) . 

IF(C0DE .NE. ' ') GO TO 500 
NPTS=NPTS+1 
100 CONTINUE 

IF( NPTS .GT. MPTS1 ) GO TO 9000 

500 CALL PCHRAP( NPTS, FLAT, FLONG,BXY) 

CALL PCAREA( BXY , NPTS , 2 , AREA, ISW) 

IF( ISW .EQ. -1) GO TO 9000 

DO 1000 1=1, NPTS 

BDCHAN( LOC 1+5+2 *I)=BXY( 1 , I) 
BDCHAN(LOCl+6+2*I)=BXY(2 , I) 

IF( INPLST .EQ. ’N’) GO TO 1000 

WRITE( IOUNIT, 2) FLAT( I ) , FLON G( I ) , BXY( 1,1), BXY( 2,1) 

1000 CONTINUE 

BDCHAN( LOC 1+4 ) =AREA 

BDCHAN( LOC 1+3 )=FLOAT( NPTS )+0.1 

CALL PCBDPK 

ALLOCATE 4-POINT TYPE 1 BOUNDARY 
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C 

MPTS2=4 

CALL PCBALC( 1 ,MPTS2 ,L0C2) 

C 

RETURN 
ERROR 

9000 WRITE( 1 ,9100) MPTS 

9100 FORMAT(/ IX, 'ERROR: THERE IS A PROBLEM WITH THE BASIN BOUNDARY', 
* /IX,' it OF BOUNDARY PTS SHOULD BETWEEN 3 AND ',13) 

C 

RETURN 1 
C 

1 F0RMAT(A1 , 2F8. 2 ) 

2 FORMAT(1X,2F10.3,10X, '( ' .2F12.3 , ' ) ' ) 

C 

END 
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C *** SUBROUTINE PCINC *** 

C 

C OBJECTIVES: 

C CONTROL GROUP 
C 1. I/O OPTIONS 

C 2. RUN TITLES 

C 3. DEBUG FLAGS 

C 
C 

SUBROUTINE PCINC( CODE , TITLE , INPLST) 

C 

CHARACTER TITLE*79 , CODE* 1 , INPLST* 1 
C 

READ( 29,1) CODE, TITLE 

I/O OPTIONS 

WRITE( 1,5) 

READ( 1,2) INPLST 
C 

IF(CODE .NE. ' ') RETURN 
READ( 29 ,2) CODE 
RETURN 
C 

C I/O OPTIONS 
C DEBUG FLAGS 
C 

1 FORMAT(Al ,A79) 

2 FORMAT(Al) 

5 F0RMAT(/1X, 'WANT INPUT DATA TO BE PRINTED? ENTER Y OR N' ) 
C 

END 
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C *** SUBROUTINE PCINN *** 

C 

C OBJECTIVES: 

C NEW DATA CASE 

C 1. CHCEK FOR VALID TYPE 0 BOUNDARY 

C 2. READ DATA AND OUTPUT OPTION 

C 3. INITIALIZE NPDTA , NLDTA , NADTA , AND NATDTA 

C 4. SET AREAL AVERAGE & COMPUTED CORRELATION OF BASIN BOUNDARY TO -1 

C 5. SET CODE INDICATING THAT DATA INPUT IS VALID, 

C I.E., PCINP, PCINL, PCINA, AND PCINI CAN BE CALLED 

C 

C 

SUBROUTINE PCINN(CODE,IOUNIT,INPLST,*) 

C 

COMMON /PCPDTA/MPDTA, NPDTA, ALPHAP( 500) ,CP( 500) ,XPDTA( 500) , 

* YPDTA( 500 ) , CAPDTA( 500 ) , VALP( 500 ) , DESCRP( 3 , 500 ) 

COMMON /PCLDTA/MLDTA, NLDTA, ALPHAL( 100) , CL( 100) , XLDTA( 2,100), 

* YLDTA( 2,100), CALDTA( 100 ) ,VALL( 100 ) , DESCRL( 3,100) 
COMMON /PCADTA/MADTA, NADTA, MATDTA, NATDTA, CA( 20) ,VALA(20) , 

* DESCRA(4 , 20) 

COMMON /PCBDCH/MCHAN , NCHAN , BDCHAN ( 1 0000 ) 

C 

CHARACTER CODE* 1 , INPLST* 1 

CHECK FOR VALID TYPE 0 BOUNDARY 

IF( IFIX(BDCHAN( 1 ) ) .NE. 0) GO TO 9000 
IF( IFIX( BDCHAN(4 ) ) .EQ. 0) GO TO 9100 

READ DATA AND OUTPUT OPTIONS 


INITIALIZE 

NPDTA=0 

NLDTA=0 

NATDTA=0 

SET AREAL AVERAGE AND COMPUTED CORRELATION OF BASIN BOUNDARY TO -1.0 

BDCHAN (.6) =-1.0 . - - _ . 

BDCHAN( 7)=-l .0 

SET CODE INDICATING THAT DATA INPUT IS VALID 


READ(29,2) CODE 
2 FORMAT(Al) 


1000 IF(CODE .EQ. 'P') 

GO 

TO 

1100 

IF(CODE .EQ. 'L') 

GO 

TO 

1200 

IF(CODE .EQ. ' A' ) 

GO 

TO 

1300 

IF(CODE .EQ. 'I') 

GO 

TO 

1400 

RETURN 
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1100 WRITE( IOUNIT, 1 1 10) 

1110 F0RMAT(/1X, 'POINT SAMPLING DATA IS READ:’) 

CALL PCINP(CODE, IOUNIT, INPLST) 

WRITE( IOUNIT, 1120) 

1120 FORMAT( IX, 'PROCESSING POINT DATA COMPLETED’) 

GO TO 1000 
C 

1200 WRITE( IOUNIT, 1210) 

1210 F0RMAT(/1X, 'LINE SAMPLING DATA IS READ:’) 

CALL PCINL(CODE, IOUNIT, INPLST) 

WRITE(IOUNIT, 1220) 

1220 FORMAT( IX, ’PROCESSING LINE DATA COMPLETED’) 

GO TO 1000 
C 

1300 WRITE( IOUNIT, 1310) 

1310 FORMAT(/ IX, ’AREAL SAMPLING DATA IS READ:’) 

CALL PCINA( CODE , IOUNIT , INPLST , * 9 200 ) 

WRITE( IOUNIT, 1320) 

1320 FORMAT( IX, ’PROCESSING AREAL DATA COMPLETED’) 

GO TO 1000 
C 

1400 WRITE( IOUNIT, 1410) 

1410 FORMAT(/ IX, ’IMAGERY SAMPLING DATA IS READ:’) 

CALL PCINI(CODE, IOUNIT, INPLST, *9200) 

WRITE( IOUNIT, 1420) 

1420 FORMAT( IX, ’PROCESSING IMAGERY DATA COMPLETED’) 

GO TO 1000 
C 

9000 WRITE( IOUNIT, 90 10) 

9010 FORMAT(/ IX, ’ERROR: BASIN BOUNDARY IS NOT OF TYPE O’) 

RE TURN 1 

9100 WRITE( IOUNIT, 91 10) 

9110 FORMAT(/ IX, ’ERROR: THERE IS NO BASIN BOUNDARY POINTS') 
RETURN1 

9200 CONTINUE 
RETURN1 
C 

END 
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C *** SUBROUTINE PCINP *** 

C 

C OBJECTIVES: 

C POINT DATA GROUP 

C 1. READ POINT DATA INTO COMMON BLOCK PCPDTA 

C 2. CONVERT (LAT. ,LONG. ) TO (X,Y) HRAP COORDINATE SYSTEMS 

C 
C 

SUBROUTINE PCINP(CODE, IOUNIT, INPLST) 

C 

COMMON /PCPDTA/MPDTA , NPDTA , ALPHAP( 500 ) , CP( 500 ) , XPDTA( 500 ) , 

* YPDTA( 500 ) , CAPDTA( 500 ) , VALP( 500 ) , DESCRP( 3,500) 

C 

DIMENSION PXY(2 ,500) 


CHARACTER CODE*l , INPLST* 1 


DO 100 1=1 .MPDTA+l 

READ( 29 , 1 ) CODE , ( DESCRP( J , I ) , J= 1 , 3 ) , XPDTA( I ) , YPDTA( I ) , VALP( I ) , 
* CP( I) , ALPHAP( I) 

IF( CODE .NE. ' ') GO TO 500 

IF( I .GT. 1 .AND. CP(I) .LT. .00001) CP( I)=CP(I-1 ) 

IF( I .GT. 1 .AND. ALPHAP(I) .LT. .00001) ALPHAP( I)=ALPHAP( 1-1 ) 
NPDTA=NPDTA+1 
100 CONTINUE 

IF( NPDTA .GT. MPDTA) GO TO 9000 


500 CALL PCHRAP ( NPDTA , XPDTA , YPDTA , PXY ) 

DO 600 1=1, NPDTA 

IF( INPLST .EQ. ' N ' ) GO TO 550 

WRITE ( IOUNIT , 2 ) (DESCRP( J,I) ,J=1 ,3) , XPDTA( I ) , YPDTA (.1 ) , VALP( I ) , 

* CP( I) , ALPHAP( I) ,PXY( 1,1), PXY(2 , I) 

550 XPDTA(I)=PXY( 1 , I) 

YPDTA( I)=PXY(2 , I) 

600 CONTINUE 

RETURN 

ERROR : -TOO MANY POINT-SAMPLES - 

9000 WRITE( IOUNIT, 9 100) NPDTA 

9100 FORMATC/1X, ’ERROR: # OF POINT-AMPLES CAN NOT EXCEED', 14, 

* /IX, 'THE LEFT OVER POINT-SAMPLES ARE IGNORED') 

C 

GO TO 500 
C 

1 F0RMAT(A1 , 3A4 , 3F8. 2 , 2F7 .4 ) 

2 F0RMAT( 1X.3A4 ,2F9.3,F11.3,2F9.4, '( ',2F9.3, ')') 

END 
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C *** SUBROUTINE PCINL *** 

C 

C OBJECTIVES: 

C LINE DATA GROUP 

C 1. READ FLIGHT LINE DATA INTO COMMON BLOCK PCLDTA 

C 2. CONVERT (LAT. ,LONG.) TO (X,Y) HRAP COORDINATE SYSTEMS 
C 
C 

SUBROUTINE PCINL(CODE , IOUNIT , INPLST) 

C 

COMMON /PCLDTA/MLDTA, NLDTA,ALPHAL( 100) ,CL( 100) ,XLDTA(2 , 100) , 

* YLDTA( 2,100), CALDTA( 100 ) ,VALL( 100 ) , DESCRL( 3 ,100) 


DIMENSION TEMPX( 500 ) , TEMPY( 500 ) , PXY( 2,500) 

C . 

CHARACTER CODE* 1 , INPLST* 1 
C 

DO 100 1= 1 , MLDTA+ 1 

READ( 29 , 1 ) CODE , ( DESCRL( J ,1) , J=1 ,3) , XLDTA( 1,1), YLDTA( 1,1), 

* XLDTA( 2,1), YLDTA( 2,1) , V ALL( I ) , CL ( I ) , ALPHAL ( I ) 
IF(CODE .NE. ’ ') GO TO 150 

IF(I .GT. 1 .AND. CL( I) .LT. .00001) CL( I)=CL( 1-1 ) 

IF( I .GT. 1 .AND. ALPHAL(I) .LT. .00001) ALPHAL(I)=ALPHAL(I-1 ) 
NLDTA=NLDTA+1 
100 CONTINUE 

IF(NLDTA .GT. MLDTA) GO TO 9000 
C 

150 IF( INPLST .EQ. 'N') GO TO 500 
DO 160 1=1 , NLDTA 

WRITE ( IOUNI T , 2 ) (DESCRL( J , I) , J=1 ,3) , XLDTA( 1,1), YLDTA( 1,1), 

* XLDTA(2 , I) ,YLDTA( 2,1) ,VALL( I) ,CL( I) , ALPHAL( I) 

160 CONTINUE 

C 

500 CONTINUE 

DO 300 J=1 ,2 

DO 200 1=1, NLDTA 
TEMPX(I)=XLDTA(J, I) 

TEMPY( I)=YLDTA( J , I) 

200 CONTINUE 
C 

CALL PCHRAP ( NLDTA , TEMPX , TEMPY , PXY) 

C 

DO 250 1=1, NLDTA 

XLDTA( J , I)=PXY( 1,1) 

YLDTA( J , I)=PXY( 2 , I) 

250 CONTINUE 
C 

300 CONTINUE 
C 

RETURN 

C 

9000 WRITE( IOUNIT, 9 100) MLDTA 

9100 FORMAT(/ IX, 'ERROR: # OF LIN-SAMPLES CAN NOT EXCEED ',14, 

* /IX,' THE LEFT OVER LINE-SAMPLES ARE IGNORED') 
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RETURN 


1 F0RMAT(A1 ,3A4,5F8.2,2F7.4) 

2 FORMAT(lX,3A4,4F9.3,F12.4,2F9.4) 
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C *** SUBROUTINE PCINA *** 

C 

C OBJECTIVES: 

C AREAL DATA GROUP 

C 1. PROCESSES ALL AREAL SAMPLE VALUES FOR A SINGLE TECHNOLOGY TO 

C PRODUCE A BASIN-WIDE AREAL AVERAGE FOR THAT TECHNOLOGY 

C 2. LOAD THE RESULTS INTO ONE ENTRY IN COMMON BLOCK PCADTA 
C 
C 

SUBROUTINE PCINA(CODE , IOUNIT , INPLST ,*) 

C 

COMMON /PCADTA/MADTA , NADTA , MATDTA , NATDTA , CA( 20 ) , VALA( 20 ) , 

* DESCRA(4 , 20) 

COMMON /PCBDCH/MCHAN ,NCHAN , BDCHAN( 10000) 

C 

DIMENSION XADTA( 4 ) , YADTA( 4 ) 

DIMENSION XYGRID(2 ,4 ) 

C 

CHARACTER CODE*! , INPLST* 1 
C 

100 CONTINUE 
C 

NADTA=0 

NATDTA=NATDTA+1 

IF (NATDTA .GT. MATDTA) GO TO 9000 
CMSUM=0. 

CDSUM=0. 

CSUM=0. 

C 

READ( 29,1) CODE , ( DESCRA( J , NATDTA) , J= 1 , 4 ) 

IF( INPLST .EQ. 'N') GO TO 300 
WRITE ( IOUNIT , 3 ) ( DESCRA( J , NATDTA) , J= 1 , 4 ) 

300 READ( 29,2) CODE , (XADTA( I) , YADTA( I) , 1=1 ,4) ,VALA1 ,TEMCA1 
IF( CODE .NE. ' ') GO TO 6000 
NADT A= NADT A+ 1 

IF (NADTA .GT. 1 .AND. TEMCA1 .LT. .00001) TEMCA1=CA1 
CA1=TEMCA1 

IF (NADTA .GT. MADTA) GO TO 9000 

CONVERT (LAT..LONG.) TO HRAP (X,Y) 

CALL PCHRAP(4 ,XADTA, YADTA,XYGRID) 

FIND SIZE OF SAMPLE AREA 

CALL PCAREA(XYGRID,4 , 2 , DAREA, ISW,*9000) 

C IF( ISW .EQ. -1) GO TO 9000 

C 

C LOAD THE AREAL SAMPLE BOUNDARY DATA INTO ARRAY BDCHAN 
C 

LAS=IFIX(BDCHAN(2)) 

C 

DO 400 J=1 ,4 

BDCHAN (LAS+5+2*J)=XYGRID( 1 , J) 

BDCHAN(LAS+6+2*J)=XYGRID( 2 , J) 

400 CONTINUE 
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C 

IF( INPLST .EQ. 'N') GO TO 370 

WRITE( I0UNIT,4) (XADTA( I) ,YADTA( I) , 1=1 ,4 ) ,VALA1 ,CA1 
WRITE(I0UNIT,5) (BDCHAN(LAS+5+2*I ) , BDCHAN(LAS+6+2*I) , 1=1 ,4 ) 
WRITE( I0UNIT,6) DAREA 
C 

370 BDCHAN ( LAS+3 ) =4 . +0 . 1 
C 

CALL PCBDPK 
C 

C FIND AUGMENTED BOUNDARIES 
C 

C 1. ALLOCATE 1ST BOUNDARY (AUGMENTED BASIN BOUNDARY) 

C 

MPTS 1 = (MCHAN-NCHAN- 1 4 ) / 6 
CALL PCBALC(2,MPTS1,L0CA1) 

2. ALLOCATE 2ND BOUNDARY (AUGMENTED AREAL SAMPLE BOUNDARY) 


MPTS2=-1 

CALL PCBALC( 2 ,MPTS2 ,LOCA2) 

3. AUGMENTED BOUNDARIES 

NF1=IFIX(BDCHAN( 4 ) ) /* # OF VERTEX IN 1ST BOUNDARY 

/ 

NF2=IFIX(BDCHAN( LAS+3)) 

CALL PCBAUG(BDCHAN(8) ,NF1 , BDCHAN (LAS+7 ) ,NF2,BDCHAN(LOCAl+7) , 

* MPTS 1 , NPTS 1 , BDCHAN(LOCA2+7 ) ,MPTS2 , NPTS2 , 

* NINTER, IER,*9100) 

BDCHAN(LOCAl+3 )=FLOAT( NPTS1 )+0. 1 
BDCHAN( LOCA2+3 ) =FLOAT( NPTS2 )+0 . 1 

IF( IER .NE. 0) GO TO 9000 
IF( NINTER .GT. 0) GO TO 2000 

NO INTERSECTION, TEST FOR WHOLE SAMPLE INSIDE BASIN 

CALL PCIFIN(BDCHAN(LAS+7 ) , BDCHAN( LAS+8 ) , BDCHAN( 8 ) , NF 1 , 2 , INOUT) 

IF( INOUT .EQ. 0) GO TO 1200 
IF( INPLST .EQ. ' N ' ) GO TO 1100 
WRITE( IOUNIT, 1050) 

1050 FORMAT( IX, 'WHOLE SAMPLE INSIDE BASIN') 

1100 CMSUM=VALA1*CA1+CMSUM 
CSUM=CA1+CMSUM 
CDSUM=CA1*DAREA+CDSUM 
GO TO 5000 

TEST FOR WHOLE BASIN INSIDE AREAL SAMPLE 
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1200 CALL PCIFIN(BDCHAN(8) , BDCHAN(9) , BDCHAN(LAS+7) , NF2 , 2 , INOUT) 
IF ( INOUT .EQ. 0) GO TO 1270 
IF( INPLST .EQ. 'N') GO TO 1250 
WRITE(IOUNIT, 1230) 

1230 F0RMAT( 1-X, 'WHOLE BASIN INSIDE SAMPLE’) 

1250 CADJ=CA1*BDCHAN(5)/DAREA 
CMSUM=CADJ*VALA1+CMSUM 
C SUM=CAD J +C SUM 
CDSUM=CADJ*DAREA+CDSUM 
GO TO 5000 

1270 IF( INPLST .EQ. ’N’) GO TO 5000 
WRITE( IOUNIT, 1280) 

1280 FORMAT( IX, 'BASIN AND AREAL SAMPLE ARE DISJOINT') 

GO TO 5000 

2000 CALL PCBDPK. 

LOCA2=IFIX(BDCHAN(LOCA1+1 ) ) 

ALLOCATE SPACE FOR SUB-BASIN 

MSUB=-1 

CALL PCBALC( 1 ,MSUB,LOCSUB) 

IF( INPLST .EQ. 'N') GO TO 2090 
WRITE ( IOUNIT, 9) 

WRITE( IOUNIT, 8) (BDCHAN(I) , 1=1 , IFIX(BDCHAN( 2 ) )-l ) 

WRITEC IOUNIT, 8) (BDCHAN(I) , I=IFIX(BDCHAN(2) ) ,L0CA1-1) 
WRITE( IOUNIT, 8) (BDCHAN( I) , I=L0CA1 , LOCA2-1 ) 

WRITE( IOUNIT, 8) ( BDCHAN ( I ) , I=LOCA2 , LOCSUB- 1 ) 


FIND OVERLAPPING AREAS 

2090 BID=0. 

ILOOP1=0 

ILOOP2=0 

2100 CALL PCSUBB(BDCHAN(L0CA1+7),NPTS1,BDCHAN(L0CA2+7),NPTS2, 

* BDCHAN ( LOCSUB+7 ) , NSUB ,MSUB , 1 , IER , 0 , *9200 ) 

IF( NSUB .EQ. 0) GO TO 2200 

BDCHAN( LOCSUB+3 ) =FLOAT( NSUB )+0. 1 
ILOOPl=l 

CALL PCAREA( BDCHAN( LOCSUB+7 ) , NSUB , 2 , AREA , I ER , *9000 ) 

BID=BID+AREA 

GO TO 2100 

2200 CALL PCSUBB(BDCHAN(LOCA2+7) ,NPTS2 ,BDCHAN(LOCAl+7) ,NPTS1 , 

* BDCHAN ( LOCSUB+7 ) , NSUB ,MSUB , 1 , IER , 0 , *9200 ) 

IF( NSUB .EQ. 0) GO TO 2300 

BDCHAN ( LOCSUB+3 ) =FLOAT( NSUB ) +0 . 1 
ILOOP2=l 

CALL PCAREA( BDCHAN( LOCSUB+7 ) , NSUB , 2 , AREA , IER , *9000 ) 

BID=BID+AREA 

GO TO 2200 
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2300 IF( IL00P1 .EQ. 1 .OR. IL00P2 .EQ. 1 ) GO TO 2500 

CALL PCSUBB( BDCHAN( L0CA1+7 ) , NPTS1 , BDCHAN(LOCA2+7 ) , NPTS2 , 

* BDCHAN(L0CSUB+7) ,NSUB,MSUB, 1 ,IER, 1 ,*9200) 

IF ( NSUB .EQ. 0) GO TO 2350 
BDCHAN( LOCSUB+3 ) =FL0AT( NSUB )+0 . 1 

CALL PCAREA( BDCHAN( LOCSUB+7 ) , NSUB , 2 , AREA , IER , *9000 ) 

BID=BID+AREA 
GO TO 2500 

2350 WRITE ( IOUNIT,2360) 

2360 F0RMAT(/1X, '* INTERSECTION PTS FOUND, BUT SUB-BASIN NOT FOUND *' 
GO TO 5000 

BID IS OVERLAPPED AREA, UPDATE SUMS 

2500 CAD J =CA1*BID/ DARE A 

CMSUM=CADJ*VALA1+CMSUM 
CSUM=CADJ+CSUM 
CDSUM=CADJ*DAREA+CDSUM 
GO TO 5000 

DEFINE VALUE FOR THIS TECHNOLOGY 

6000 IF(CSUM .EQ. 0) GO TO 6100 
VALA( NATDTA) =CMSUM/ CSUM 
CA( NATDTA) =CDSUM/BDCHAN( 5 ) 

GO TO 5000 

6100 VALA( NATDTA) =0. 

CA( NATDTA) =CDSUM/ BDCHAN( 5 ) 

- GO TO 5000 

CLEANUP: JUST DELETES SPACE ALLOCATED FOR AUGMENTED BOUNDARIES 

5000 IF(L0CA1 .GT. 0) BDCHAN( LOCA 1 ) =- 1 . 0 
IF(LOCA2 .GT. 0) BDCHAN(L0CA2)=-1 .0 
IF(L0CSUB .GT. 0) BDCHAN(LOCSUB)=-l .0 

CALL PCBDPK 


IF(C0DE .EQ. ' ') GO TO 300 
IF (CODE .EQ. ’A’) GO TO 100 
RETURN 

ERROR 
C 

9000 WRITE(IOUNIT, 9010) . 

9010 FORMAT(/ IX,. 'THERE IS A PROBLEM CALLING PCAREA IN AREAL SAMPLES') 
RETURN1 
C 

9100 WRITE ( IOUNIT.91 10) 

9110 F0RMAT(7 IX, ’ THERE IS A PROBLEM -IN CALLING-PCBAUG:' ); 

RETURN1 
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9200 WRITE(I0UNIT,9210) 

9210 FORMAT(/ IX, 'THERE IS A PROBLEM IN CALLING PCSUBB') 
RETURN1 
C 

1 F0RMAT(A1,4A4) 

2 F0RMAT(A1,9F8.2,F7.4) 

3 F0RMAT(/1X,4A4) 

4 F0RMAT(/1X,9F8.2,F7.4) 

5 FORMAT( ' ( ' ,8F8. 2 , ' ) ' ) 

6 F0RMAT(1X, 'AREA = ',F14.4) 

8 F0RMAT(1X,8F9.2) 

9 FORMAT( IX, 'AUGMENTED BOUNDARIES:') 

C 

END 
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C *** SUBROUTINE PCINI *** 

C 

C OBJECTIVES: 

C IMAGE DATA GROUP 

C 1. PROCESSES ALL IMAGERY SAMPLE VALUES FOR A SINGLE TECHNOLOGY TO 
C PRODUCE AN INAGERY AVERAGE FOR THAT TECHNOLOGY 

C 2. LOADS THE RESULTS INTO ONE ENTRY IN COMMON BLOCK PCADTA 
C 
C 

SUBROUTINE PCINI(CODE , IOUNIT , INPLST,*) 

C 

COMMON /PCADTA/MADTA , NADTA , MATDTA , NATDTA , CA(20 ) , VALA( 20 ) , 

* DESCRA(4 , 20) 

COMMON /PCBDCH/MCHAN , NCHAN ,BDCHAN( 10000) 

C 

DIMENSION PV(64 ) 

CHARACTER* 1 PVCC( 64 ), CHARC( 64 ), INPLST* 1 
DIMENSION XYGRID( 2 , 1 ) 

C 

CHARACTER CODE*l 
C 

100 CONTINUE 
NADTA=0 

NATDTA= NATDTA+ 1 

IF (NATDTA .GT. MATDTA) GO TO 9000 
C 

READ( 29,1) CODE , ( DESCRA( J , NATDTA) , J= 1 , 4 ) , C 
READ( 29,3) CODE , NOVAL , ( PVCC( I ) , PV ( I ) , 1= 1 , NOVAL) 

C 

IF( INPLST .EQ. 'N') GO TO 200 

WRITE ( IOUNIT , 2 ) ( DESCRA( J , NATDTA) , J= 1 , 4 ) , C 

WRITE( IOUNIT , 4 ) NOVAL , ( PVCC( I ) , PV( I ) , 1=1 , NOVAL) 

C 

200 SUMPV=0. 

NPIB=0 /* 0 OF PIXELS IN BASIN */ 

C 

300 READ(29 ,6) CODE ,XST, YST, XEND, YEND, NPL, (CHARC( I) , 1=1 , NPL) 

C 

IF( CODE .NE. ' ') GO TO 6000 
NADT A= N ADTA+ 1 

IF( NADTA .GT. MADTA) GO TO 9000 
C 

IF( INPLST .EQ. ’N') GO TO 400 

WRITE ( IOUNIT , 7 ) XST , YST , XEND , YEND , NPL , ( CHARC( I ) , 1= 1 , NPL) 

CONVERT (LAT. ,LONG. ) TO HRAP (X,Y) 

400 CALL PCHRAP( 1 ,XST, YST.XYGRID) 

XST=XYGRID( 1,1) 

. YST=XYGRID( 2,1) 

C 

CALL PCHRAP(1 , XEND, YEND, XYGRID) 

XEND=XYGRID( 1,1) 

YEND=XYGRID( 2,1) 
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C 

IF( INPLST .EQ. 'N') GO TO 450 
WRITE(I0UNIT,8) XST, YST.XEND, YEND 
C 

450 DX=(XEND-XST)/(NPL-1) 

DY=( YEND— YST) / ( NPL-1 ) 

PROCESS A LINE OF IMAGE 

DO 500 1=1, NPL 
XI=XST+DX* ( I- 1 ) 

YI=YST+DY*( 1-1 ) 

DO 550 J=i , NOVAL 

IF(CHARC( I) .EQ. PVCC(J)) GO TO 560 
550 CONTINUE 

GO TO 9000 
560 PVI=PV(J) 

NF1=IFIX( BDCHAN(4 ) ) 

CALL PCIFIN(XI , YI , BDCHAN( 8) ,-NFl , 2 , INOUT) 

IF( INOUT .LT. 1) GO TO 500 
SUMPV=SUMPV+PVI 
NPIB=NPIB+1 
500 CONTINUE 

GO TO 300 

DEFINE VALUE FOR THIS TECHNOLOGY 

6000 IF( NPIB .EQ. 0) GO TO 6100 
VALA( NATDTA)=SUMPV /NPIB 
CA(NATDTA)=C 
C 

IF(CODE .EQ. 'I’) GO TO 100 /* ANOTHER TECHNOLOGY */ 

RETURN 
C 

6100 VALA( NATDTA) =0 . 

CA(NATDTA)=0. 

IF(CODE .EQ. 'I') GO TO 100 
RETURN 
C 

1 F0RMAT(A1,4A4,F7.4) 

2 F0RMAT(/1X,4A4,F7.4) 

3 F0RMAT(A1 , 14 ,8(A1 ,F8. 2) ,7(/5X,8(Al ,F8. 2)) ) 

4 F0RMAT(1X,I4 ,8(A1 ,F8. 2) ,7 (/5X, 8(A1 ,F8 . 2) ) ) 

-6 F0RMAT(A1 ,4F8. 2 , 14 ,43A1 ) 

7 F0RMAT(1X,4F8.2,I4,43A1) 

8 FORMAT( '( ',4F8.2, ') ') 

9000 WRITE(IOUNIT,9100) 

9100 FORMAT(/ IX, 'THERE IS A PROBLEM IN IMAGERY SAMPLES') 

C 

RETURN 1 
END 
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*** SUBROUTINE PCSUBB *** 

OBJECTIVES: 

PROCESSES TWO AUGMENTED BOUNDARY SETS (FROM SUBROUTINE PCBAUG) TO 
FIND ONE SUB-BASIN 

TO FIND ALL SUB-BASINS INSIDE BOTH BOUNDARIES, PCSUBB MUST BE 
CALLED REPEATEDLY UNTIL NSUBB=0 

PRINCIPAL VARIABLES: 


NAME 

I/O 

DIMENSION 

DESCRIPTION 

G1 

I 

(3.NG1) 

AUGMENTED BOUNDARY 1 

NG1 

I 

- 

if OF VERTICES IN G1 

G2 

I 

(3 , NG2) 

AUGMENTED BOUNDARY 2 

NG2 

I 

- 

if OF VERTICES IN G2 

SUBB 

0 

(2, MSUBB) 

VERTEX PTS OF A SUB-BASIN INSIDE BOTH BOUNDARIES 
1 AND 2 (FOR INOUT=l ) 

FOR IN0UT=0, SUBB DEFINES A SUB-BASIN INSIDE 
BOUNDARY 1 BUT NOT INSIDE BOUNDARY 2 

NSUBB 

0 

- 

if OF VERTICES IN SUBB 

MSUBB 

I 

- 

MAX if OF VEPTHCES IN SUBB 

INOUT 

I 


=1 (NORMAL) FIND A SUB-BASIN INSIDE BOTH 1 & 2 
=0 TO FIND A SUB-BASIN INSIDE 1 BUT OUTSIDE 2 

IER 

0 

— 

ERROR FLAG , =0 NORMAL RETURN 

=1 ERROR - NOT ENOUGH SPACE 

IVPIB 

I 

- 

=0 SOME VERTEX-PTS INSIDE BOUNDARY OF OTHER 


=1 NO VERTEX-PT -INSIDE BOUNDARY OF EACH OTHER 


SUBROUTINES CALLED: 
PCIFIN 


SUBROUTINE PCSUBB(G1 , NG1 ,G2 , NG2 , SUBB , NSUBB ,MSUBB , INOUT, IER, 

* IVPIB,*) 

DIMENS ION G1 ( 3 , NG1 ) , G2 ( 3 , NG2 ) , SUBB( 2 ,MSUBB) 

IER=0 
NSUBB=0 

FIND A STARTING POINT 

IF( IVPIB .EQ. 1) GO TO 511 
DO 500 IST=1 , NG1 

IF(ABS(G1 ( 3 , 1ST) ) .GT. 0.1) GO TO 500 /* INTERESCTED OR USED 

*/ 

CALL PCIFIN(G1(1,IST),G1(2,IST),G2, NG2 ,3 , ITST) 

IF( ITST .EQ. INOUT) GO TO 550 
500 CONTINUE 
C 

RETURN 
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511 DO 521 IST=1,NG1 
IEND=IST+1 

IF( IEND .GT. NG1 ) IEND=1 

IF(G1 ( 3 , 1ST) .GT. 0.1 .AND. G1(3,IEND) .GT. 0.1) GO TO 531 
521 CONTINUE 

RETURN 

531 ISTG2=IFIX(G1 (3 , 1ST) ) 

START POINT FOUND 

550 I=IST 
ISET=1 

ADD POINT I 

600 IF(NSUBB .EQ. MSUBB) GO TO 9000 
NSUBB=NSUBB+1 
DO 1000 J=1 ,2 

IF( ISET .EQ. 1) SUBB( J ,NSUBB)=G1 ( J , I) 

IF( ISET .EQ. 2) SUBB(J,NSUBB)=G2(J, I) 

1000 CONTINUE 

NOW MARK THIS POINT AS USED 

IF( ISET . NE. 1) GO TO 1050 
IF(G1 (3 , I) .LT. 0.1) G 1 ( 3 , I ) — — 1.0 
GO TO 1100 

1050 IF( ISET .NE. 2) GO TO 1150 

IF(G2(3, I) .LT. 0.1) G2 ( 3 , I ) =- 1 . 0 

1100 IF( ISET .EQ. 2) GO TO 1200 

FOR SET 1 

1150 1=1+1 

IF( I .GT. NG1 ) 1=1 
IF( I .EQ. 1ST) RETURN 

CHECK FOR INTERSECTION 

IF(G1(3,I) .LT. 0.1) GO TO 600 
ISET=2 

I=IFIX(G1(3 , 1) ) 

GO TO 600 

FOR SET 2 

1200 1=1+1 

IF( INOUT .EQ. 0) 1=1-2 
IF(I .LT. 1) I=NG2 
IF( I .GT. NG2 ) 1=1 

IF( IVPIB .EQ. 1 .AND. I .EQ. ISTG2 ) RETURN 
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CHECK FOR INTERSECTION 
C 

IF(G2(3,I) .LT. 6.1) GO TO 600 
ISET=1 

I=IFIX(G2(3,I)) 

GO TO 600 
C 
C 

C ERPOR : NOT ENOUGH SPACE FOR SUB-BASIN POINT 
C 

9000 IER=-1 
RETURN 1 
END 



non onnoo 


C *** PCBAUG *** 

C 

C OBJECTIVES: 

C PRODUCES THE AUGMENTED BOUNDARY SETS WHICH INCLUDE ALL INTERSECTION 
C POINTS OF TWO BOUNDARIES 
C 

C PRINCIPAL VARIABLES: 


c 

NAME 

I/O 

DIMENSION 

DESCRIPTION 

L# 

c 

F1XY 

I 

( 2 , NF 1 ) 

X-Y PAIRS OF FIRST BOUNDARY (1=X, 2=Y) 

c 

NF1 

I 

- 

it OF PARIR IN FIRST BOUNDARY 

c 

F2XY 

I 

(2 , NF2) 

X-Y PAIRS OF SECOND BOUNDARY 

c 

NF2 

I 

- 

It OF PAIRS IN SECOND BOUNDARY 

c 

G1XY 

0 

(3 ,MG1) 

AUGMENTED FIRST BOUNDARY INCLUDING INTERSECTIONS 

c 




1=X,2=Y,3=FLAY 

c 

MG1 

I 

- 

MAX. it OF POINTS IN AUGMENTED FIRST BOUNDARY 

c 

NG1 

0 

- 

ACTUAL it OF PTS IN COMPLETED FIRST BOUNDARY 

c 

G2XY 

0 

(3 ,MG2) 

SECOND AUGMENTED BOUNDARY 

c 

MG2 

I 

- 

MAX. it OG POINTS 

c 

NG2 

0 

- 

ACTUAL it OF PTS 

c 

NINTER 

0 

- 

it OF INTERSECTIONS BETWEEN F1XY AND F2XY 

c 

IER 

0 

- 

RETURN STATUS, =0 NO ERROR 

c 




=1 ERROR OCCURRED (TOO LITTLE 

STORAGE) 





SUBROUTINE PCBAUG( F 1XY , NF1 , F2XY , NF2 , G1 XY ,MG1 , NG1 , G2XY ,MG2 , NG2 , 

NINTER, IER ,*) 


DIMENSION F1XY( 2 ,NF1 ) , F2XY( 2 ,NF2) , G1XY( 3 ,MG1 ) , G2XY( 3 ,MG2) 
DIMENSION G2XY0( 1000) 


IER=0 
NINTER=0 
DO 200 1=1, NF1 
DO 100 J=1 ,2 

100 G1XY(J,I)=F1XY(J,I) 

200 G1XY(3,I)=0. 

DO 250 1=1 , NF2 
DO 150 J=1 ,2 

150 G2XY( J , I)=F2XY( J , I) 

250 G2XY( 3 , I)=0. 

NG1=NF1 

NG2=NF2 

FIND ALL INTERSECTIONS, ADD TO G1XY & G2XY 
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C 

1000 


C 


1200 

C 

1300 


C 

ADD 
C 

1500 

C 

1600 


CONTINUE 

DO 1900 IS=1 ,NG1 

IE=IS+1 /* ENDING POINT OF A LINE */ 

IF(IE .GT. NG1 ) IE=1 
X1S=G1XY( 1 , IS) 

X1E=G1XY( 1 , IE) 

IFL1S=IFIX(G1XY( 3 , IS) ) 

Y1S=G1XY(2 , IS) 

Y1E=G1XY( 2 , IE) 

IFL1E=IFIX(G1XY(3,IE)) 

A1=X1S-X1E 

B1=Y1S-Y1E 

C1=X1S*Y1E-Y1S*X1E 


DO 1800 JS=1 , NG2 
J£=JS+1 

IF( JE .GT. NG2 ) JE=1 

IF( IFL1 S .EQ. 0 .AND. IFL1E .EQ. 0) GO TO 1300 
IFL2S=IFIX(G2XY( 3 , JS) ) 

IFL2E=IFIX(G2XY( 3 , JE) ) 

IF( IFL1S .EQ. 0) GO TO 1200 

IF(IFL1S .EQ. IFL2S .OR. IFL1S .EQ. IFL2E) GO TO 1800 
IF( IFL1E .EQ. 0) GO TO 1300 

IF( IFL1E .EQ. IFL2S .OR. IFL1E .EQ. IFL2E) GO TO 1800 

X2S=G2XY( 1 , JS) 

X2E=G2XY( 1 , JE) 

Y2S=G2XY( 2 , JS) 

Y2E=G2XY(2 , JE) 

A2=X2S-X2E 
B2=Y2S-Y2E 
C2=X2S*Y2E-Y2S*X2E 
D=B2*A1-B 1*A2 


IF(D .EQ. 0.0) GO TO 1800 
XSTAR=(C1*A2-C2*A1 )/D 
IF((XSTAR-X1S)*(XSTAR-X1E) .GT. 
IF( (XSTAR-X2S)*(XSTAR-X2E) .GT. 
YSTAR=(C1*B2-C2*B1 )/D 
IF( ( YSTAR-Y1S)*( YSTAR-Y1E) .GT. 
IF( ( YSTAR-Y2S)*( YSTAR-Y2E) .GT. 


0.0) GO TO 1800 
0.0) GO TO 1800 

0.0) GO TO 1800 
0.0) GO TO 1800 


POINT (XSTAR, YSTAR) TO G1XY.G2XY AND RE-ORDERING 

IF(NG1 .GE. MG1 .OR. NG2 .GE. MG2 ) GO TO 1700 

DO 1500 I=NG1 , IS ,-l 
DO 1500 J=1 ,3 

G1XY( J , 1+1 )=G1XY( J , I) 

CONTINUE 


DO 1600 I=NG2 , JS ,-l 
DO 1600 J=1 ,3 

G2XY( J , I'+l ) =G2XY( J , I ) 
CONTINUE 
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C 


NINTER=NINTER+1 
FLAG=FLOAT(NINTER)+0. 1 
G1XY( 1 , IS+1 )=XSTAR 
G1XY(2,IS+1)=YSTAR 
G1XY(3,IS+1) =FLAG 
G2XY( 1 , JS+1 )=XSTAR 
G2XY(2 , JS+1 )=YSTAR 
G2XY( 3 , JS+1 )=FLAG 
NG1=NG1+1 
NG2=NG2+1 
GO TO 1000 
C 

1700 XER=1 

RETURN1 

C 

1800 CONTINUE 
1900 CONTINUE 

RENUMBER FLAGS OF ALL INTERSECTION POINTS 

IF(NINTER .EQ. 0) RETURN 

DO 2000 K=1 , NG2 
2000 G2XY0(K)=G2XY(3,K) 

DO 2500 1=1, NG1 

IFL1=IFIX(G1XY(3,I)) 

IF( IFL1 .EQ. 0) GO TO 2500 
DO 2400 J-1.NG2 

IFL2=IFIX(G2XY0( J) ) 

IF( IFL2 .NE. IFL1) GO TO 2400 
G1XY(3,I)=FLOAT(J)+0. 1 
G2XY ( 3 , J ) =FLOAT ( I ) +0 . 1 

GO TO 2500 
2400 CONTINUE 


IF FULL THROUGH J LOOP: ERROR, DID NOT FIND INTERSECTION PT IN 2ND LIST 
WRITE( 1,1) 

1 FORMAT(/ IX, 'ERROR: DID NOT FIND INTERSECTION PT IN 2ND LIST') 

RETURN 1 

500 CONTINUE 


RETURN 

END 
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C *** SUBROUTINE PCINE *** 

C 

C OBJECTIVES: 

C END OF DATA CASE 

C 1. COMBINES ALL AREAL AND IMAGE TECHNOLOGIES IN COMMON BLOCK PCADTA 
C INTO A SINGLE BASIN ESTIMATE (BY CALLING PCCOMB) 

C 2. USES GRID METHOD TO COMBINE ALL SAMPLING TECHNOLOGIES (POINT, LINE 
AND AREAL) (BY CALLING PCGRID) 

3. OUTPUTS RESULTS 

4. SETS CODE INDICATING THAT DATA INPUT IS NOT VALID UNTIL AFTER A 
CONTROL CODE 'N' 


SUBROUTINE PCINE( CODE , TITLE , IOUNIT , * ) 

C 

COMMON /PCPDTA/MPDTA , NPDTA , ALPHAP( 500 ) , CP( 500 ) , XPDTA( 500 ) , 

* YPDTA(500) , CAPDTA( 500) ,VALP(500) ,DESCRP(3 , 500) 
COMMON /PCLDTA/MLDTA , NLDTA , ALPHAL ( 100) ,CL( 100) ,XLDTA(2 , 100) , 

* YLDTA(2, 100) ,CALDTA( 100) ,VALL( 100) ,DESCRL(3, 100) 
COMMON /PCADTA/MADTA , NADTA , MATDTA , NATDTA , CA( 20 ) , VALA( 20 ) , 

* DESCRA(4 , 20) 

COMMON /PCBDCH/MCHAN , NCHAN , BDCHAN ( 10000) 

C 

CHARACTER CODE* 1 , TITLE*79 
C 

C COMBINE ALL EXISTING AREAL AND IMAGE TECHNOLOGIES IN COMMON BLOCK 
C PCADTA INTO A SINGLE BASIN ESTIMATE 
C 

CALL PCCOMB( IOUNIT) 

C 

WRITE(1,10) 

10 F0RMAT(/1X, 'ENTER // OF GRID POINTS (INTEGER IN FOUR DIGITS)') 

11 READ( 1,12) NGRID 

12 FORMAT( 14 ) 

C 

IF(NGRID .LE. 0) GO TO 1000 

CALL PCGRID( NGRID , AVAL , ACORR , TITLE , IOUNIT , *9900 ) 

C 

WRITE (IOUNIT, 5) AVAL, ACORR 

5 F0RMAT(//1X, 'THE BASIN-WIDE AREAL AVERAGE ESTIMATE = ',F12.4, 

* /IX, 'THE OVERALL ACCURACY OF THIS ESTIMATE = \F6.4) 

C - - 

WRITE( 1,14) 

14 F0RMAT(//1X, 'ENTER ANOTHER // OF GRID POINTS (IF NO MORE,' 

* ' ENTER 0)') 

GO TO 11 

C 

1000 READ( 29,1) CODE 
1 FORMAT(Al) 

IF(CODE .EQ. 'C' .OR. CODE .EQ. 'N' .OR. CODE .EQ. ’S') RETURN 
WRITE( IOUNIT, 7) 

7 F0RMAT(//1X, 'DATA CASE INPUT IS NOT VALID’) 

9900 CONTINUE 

RETURN 1 - 

C 

END 
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C ** SUBROUTINE PCCOMB *** 

C 

C OBJECTIVES: 

C COMBINES ALL AVAILABLE AREAL AND IMAGERY SAMPLING TECHNOLOGIES 
C INTO A SINGLE BASIN-WIDE AREAL ESTIMATE 
C 
C 

SUBROUTINE PCCOMB( IOUNIT) 

C 

COMMON /PCBDCH/MCHAN, NCHAN ,BDCHAN( 10000) 

COMMON / PC ADTA/MADTA , NADTA , MATDTA , NATDTA ,CA(20),VALA(20), 

* DESCRA(4 ,20) 

INITIALIZE 

BDCHAN(6)=-1 .0 
BDCHAN(7)=0.0 

IF (NATDTA .LT. 1) RETURN 
WRITE( IOUNIT, 1) 

SUM1=0.0 
SUMBM=0.0 
SUM2=0.0 

FIND A FACTOR 

DO 500 1=1, NATDTA 

IF(CA( I) ,GT. 0.999) CA(I)=.999 
IF(CA( I) .LT. 0.001) CA( I)=. 001 
SUM1=SUM1+(CA( I)/ ( l-CA(I) ) ) 

500 CONTINUE 

WEIGHT EACH MEASUREMENT 

DO 600 1=1, NATDTA 

BI=CA(I)/(SUM1*( l.-CA(I))) 

SUMBM=BI*VALA( I)+SUMBM 
SUM2=BI/ SUM1+SUM2 

WRITE (IOUNIT, 3) I ,CA( I) , VALA( I) , BI 
600 CONTINUE 

COMBINE 

BDCHAN( 7 )=SUMBM 
BDCHAN(6)=1 .0/(1 .+SUM2) 

C 

1 FORMAT(/ IX, 'AREAL SAMPLING SUMMARY:', 

* /3X, 'SAMPLING TECH CORRELATION VALUE WEIGHT', 

* /3X, ' ’) 

3 FORMAT(/8X,I2,llX,F7.4,5X,F8.2,4X,F7.4) 

RETURN 

END 
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*** SUBROUTINE PCSYMBO *** 


SUBROUTINE PCSYMBO( NX , IHTYP , IHNUM , IOUNIT) 
DIMENSION IHTYP( NX) , IHNUM( NX) 

CHARACTER*! TYPE( 130) , LCHAR( 1 30 ) * 1 , NCHAR( 130)*! 


DATA LCHAR/ , A', , B', , C', , D , ,'E', , F , , , G','H , > 'I', 

* 'NVOVPVQ'.'RVSVTVUVV'.'W, 

DATA NCHAR/ ’ 1 2 ' , '3' ,‘4’ , • 5 ' ,'6', ' 7 ' ,'8*, *9', 

* Ill*’ %'/ 


DO 9000 1=1, NX 

IF( IHTYP( I ) .NE. 0) GO TO 5000 
TYPE( I)= ' * * 

GO TO 9000 
5000 CONTINUE 

GO TO (1000,2000,3000), IHTYP(I) 
C 

1000 TYPE( I)= ' $ ' 

GO TO 9000 
C 

2000 TYPE( I)=LCHAR( IHNUM( I) ) 

GO TO 9000 
C 


3000 TYPE ( I ) =NCHAR( IHNUM( I) ) 
C 


9000 CONTINUE 


WRITE (IOUNIT, 9900) (TYPE(K) ,K=1 ,NX) 
9900 FORMAT( IX, 130A1 ) 

RETURN 

END 


•J 

'X 

'0 


'K'.'L'.'M', 
' Y ' , 105* ' Z ' / 
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C 

C 

C 

C 

C 

C 

c 

c 


c 

c 

c 


*** SUBROUTINE PCGRID *** 

OBJECTIVES: 

USING GRID METHOD, COMBINES ALL DATA (POINT, LINE, AREAL, IMAGERY) 
INTO A SINGLE BASIN-WIDE ESTIMATE 


SUBROUTINE PCGRID( NGRID , AVAL , ACORR , TITLE , IOUNIT , *) 

COMMON /PCPDTA/MPDTA , NPDTA , ALPHAP( 500 ) , CP( 500 ) ,XPDTA(500) , 

YPDTA( 500 ) , CAPDTA( 500 ) , VALP( 500 ) , DESCRP( 3,500) 
COMMON /PCLDTA/MLDTA , NLDTA , ALPHAL( 100) ,CL( 100) , XLDTA( 2 , 100) , 

YLDTA( 2,100), CALDTA( 1 00 ) ,VALL( 100), DESCRL( 3,100) 
COMMON / PCADTA/MADTA , NADTA , MATDTA , NATDTA , CA( 20 ) , VALA( 20 ) , 
DESCRA(4 , 20) 

COMMON /PCBDCH/MCHAN , NCHAN , BDCHAN( 10000 ) 

CHARACTER ANSWER*! ,TITLE*79 ,LCHAR( 130)* 1 ,NCHAR( 130)*1 


DIMENSION IHTYP( 130) , IHNUM( 130) 


DATA LCHAR/'A', 

'B' 

,'C' 

' , 'D ','E ’,’F VG’ , 'H VI' 

* ' N * , ' 0 ' , 

t p . 

,'Q' 

1 , 'R' , 'S' , 'T' , 'U' , 'V' , 'W' 

DATA NCHAR/ ' 1 ' , 

'2' 

,'3' 

, ’4’ , '5* , ’6’ , ’7' , '8' , '9' 

* ' + ’, 

f _ » 

» — 1 
y 



K\ 'L', 'M', 
Y ' , 105* 'Z ' / 

#V&', 


VALIDITY TESTS 

IF( NGRID .LT. 100) GO TO 9000 
80 IF(BDCHAN(4) .LT. 3) GO TO 9100 

IF(NPDTA+NLDTA+NATDTA .LT. 1) GO TO 9200 
C 

C DEFINE GRID 

C FIND THE WESTMOST, EASTMOST, SOUTHMOST, AND NORTHMOST POINTS OF THE 
C BASIN BOUNDARY RESPECTIVELY 
C 

XMIN=BDCHAN(8) 

XMAX=BDCHAN(8) 

YMIN=BDCHAN(9) 

YMAX=BDCHAN ( 9 ) 

DO 100 1=1 ,BDCHAN(4) 

XMIN=AMIN1 (XMIN , BDCHAN( 6+2*1) ) 

XMAX=AMAX 1 ( XMAX , BDCHAN ( 6+2 * I ) ) 

YMIN=AMIN1(YMIN,BDCHAN( 7+2*1)) 

YMAX= AMAX 1 ( YMAX , BDCHAN ( 7+ 2 * I ) ) 

100 CONTINUE 
C 

XLEN=XMAX-XMIN 
YLEN=YMAX-YMIN 
DE LTA=S QRT ( XLE N* YLE N/ NGRID) 

NX=IFIX(XLEN/DELTA+0 .5) 

NY=IFIX(YLEN/DELTA+0.5) 

NPB=BDCHAN(4) 

C 

C INITIALIZE 
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NGIB=0 

/* 

# OF GRIDS IN BASIN */ 


CSUM=0.0 

/* 

SUM OF GRIDS WTS */ 


DO 600 1=1, NLDTA 



600 

CALDTA( I)=0 .0 

/* 

FOR ALL LINE SAMPLES */ 


DO 700 1=1 , NPDTA 



700 

CAPDTA( I)=0.0 

/* 

FOR ALL POINT SAMPLES */ 


CAASAM=0.0 

/* 

SUM OF G.W. FOR AREAL TECH 


c 

C GRID POINT LOOP 
C 


WRITE( IOUNIT,800) 

800 FORMAT(// IX, 'SAMPLING TECH DISTRIBUTION IN THE BASIN',/) 
C 

DO 2900 IY=1 , NY 

YGRID=YMAX-0 . 5*DELTA-DELTA*( IY-1 ) 

DO 2800 IX=1 , NX 

XGRID=XMIN+0 . 5*DELTA+DELTA*( IX-1 ) 

C 

C TEST FOR INSIDE BASIN 
C 

IHTYP( IX)=0 

CALL PCIFIN(XGRID , YGRID , BDCHAN(8) , NPB, 2 , INOUT) 

IF( INOUT .EQ. 0) GO TO 2800 
C 

C AREAL DATA TESTING 
C 

CHIGH=BDCHAN( 6 ) . 

IHTYP(IX)=1 

IHNUM(IX)=0 

C 

C POINT DATA (CHECK FOR THE MOST HIGHLY CORRELATED SAMPLE) 

C 

DO 2100 1=1 ,NPDTA 
XOF=XGRID-XPDTA( I ) 

YOF=YGRID-YPDTA( I ) 

D=SQRT( XOF*XOF+YOF*YOF) 

C=CP( I)*EXP(-ALPHAP( I) *D) 

IF(C .LT. CHIGH) GO TO 2100 
CHIGH=C 
IHTYP( IX)=2 
IHNUM( IX)=I 
2100 CONTINUE 
C 

C LINE DATA (CHECK FOR THE MOST HIGHLY CORRELATED SAMPLE) 

C 

DO 2250 1=1 , NLDTA 
IF(CL( I) .LT. CHIGH) GO TO 2250 
X1=XLDTA( 1 ,1) 

Y1=YLDTA( 1,1) 

X2=XLDTA(2 , I) 

Y2=YLDTA( 2,1) 

R1SQ=(XGRID-X1 )**2+( YGRID-Y1 )**2 
R2SQ=(XGRID-X2)**2+(YGRID-Y2)**2 
R3SQ=(X1-X2 )**2+( Y1-Y2 )**2 
RSQMX=AMAX1(R1SQ,R2SQ) 
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RSQMN=AMIN1 (R1 SQ , R2SQ) 

IF(RSQMX .LT. (RSQMN+R3SQ) ) GO TO 2210 /* CLOSER TO LINE */ 

D=SQRT(RSQMN) /* CLOSER TO END PT */ 

GO TO 2220 

2210 R1=SQRT(R1SQ) 

R2-SQRT(R2SQ) 

R3=SQRT(R3SQ) 

S=(Rl+R2+R3)/2 

D=2.*SQRT(S*(S-R1 )*(S-R2)*(S-R3))/R3 
2220 C=CL ( I ) *EXP( -ALPHAL ( I ) *D ) /* DISTANCE FOUND */ 

IF(C .LT. CHIGH) GO TO 2250 
CHIGH=C 
IHTYP( IX)=3 
IHNUM( IX)=I 

2250 CONTINUE « 

C 

C MOST HIGHLY CORRELATED SAMPLE NOW FOUND 
C 

NGIB=NGIB+1 

CSUM=CSUM+CHIGH 

C 

GO TO (2301, 2302, 2303), IHTYP(IX) 

C 

2301 CAASAM=CAASAM+CHIGH 
GO TO 2800 

C 

2302 CAPDTA( IHNUM( IX) )=CAPDTA( IHNUM( IX) )+CHIGH 
GO TO 2800 

C 

2303 CALDTA( IHNUM( IX) )=CALDTA( IHNUM( IX) )+CHIGH 
C 

C THIS GRID POINT FINISHED 
C 

2800 CONTINUE 
C 

CALL PCSYMBO( NX , IHTYP , IHNUM , IOUNIT) 

C WRITE( IOUNIT, 2950) ( IHTYP(K) ,K=1 , NX) 

C 

2900 CONTINUE 
C 

C2950 FORMAT( IX, 130A1) 

C 

C ALL GRID POINTS FINISHED 
C 

WRITE( IOUNIT, 5556) NGIB 
WRITE( IOUNIT, 5558) DELTA 

5558 F0RMAT(1X, 'THE WIDTH OF THE GRID = ' ,F12.4) 

5556 F0RMAT(/1X, '// OF GRID POINTS ON THE BASIN = ',15) 

C 

ACORR=CSUM/NGIB 

C 

C FOR THOSE GRID PTS THAT ASSOCIATED WITH AREAL SAMPLE 
C 

AVAL=BDCHAN( 7 )*CAASAM/CSUM 
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WRITE (IOUNIT, 11) TITLE 
11 FORMAT (//1X,A79, 

* // IX, 2X, 'SAMPLING TECH' ,3X, 'MEASUREMENT VALUE' ,3X, 'WEIGHT' , 

* 3X, ' CORRELATION ' , 3X , ' ALPHA ' , 3X , ' SYMBOL ' ,. 

* /IX, 2X, ' ' , 3X, 1 ' ,3X, ' ', 

* 3X, ' ’ ,3X, ' ' ,3X, ' ',/) 

WRITE( IOUNIT ,15) BDCHAN( 7 ) , CAASAM/CSUM, BDCHAN( 6 ) 

15 FORMAT( IX, 6X, ' A' , 13X,F12 .4 , 6X,F6 . 4 , 5X,F6. 3 , 5X,6X,4X, ' ( $) ' , / ) 

16 FORMAT ( IX, 6X, ' P ' ,13X,F12.4,6X,F6.4,5X,F6.3,6X,F5.3,4X, '( ' ,A1, ') ') 

17 FORMAT( IX, 6X, 'L ' , 13X.F12 .4 ,6X ,F6 .4 , 5X,F6. 3 ,6X,F5 . 3 ,4X, ' ( ' ,A1, ') ') 


FOR THOSE GRID PTS THAT ASSOCIATED WITH POINT SAMPLE 

DO 3000 1=1 , NPDTA 

AVAL=AVAL+VALP( I)*CAPDTA( I) / CSUM 

WRITE( IOUNIT , 16)VALP( I) ,CAPDTA( I)/CSUM,CP( I) ,ALPHAP( I) , 

* LCHAR(I) 

3000 CONTINUE 

WRITE ( IOUNIT, 3010) 

3010 FORMAT( IX, ' ') 

FOR THOSE GRID PTS THAT ASSOCIATED WITH LINE SAMPLE 

DO 4000 1= 1 , NLDTA 

AVAL=AVAL+VALL( I )*CALDTA( I ) /CSUM 

WRITE( IOUNIT , 1 7 ) VALL( I ) , CALDTA( I ) / CSUM , CL( I ) , ALPHAL( I ) , 

* NCHAR(I) 

GO TO 4000 

4000 CONTINUE 

RETURN 

9000 WRITE( 1,9010) NGRID 

9010 F0RMAT(/1X, 'NO. OF GRID POINTS IS ',14,' NOT LARGE ENOUGH', 

* /IX, 'DO YOU WANT TO CONTINUE? ENTER " " IF YES') 
READ( 1,9015) ANSWER 

9015 FORMAT(Al) 

IF( ANSWER .EQ. ' ') GO TO 80 
RETURN1 


C 

C 


9100 WRITE ( IOUNIT, 9 1 10) 

9110 FORMAT(/ IX, 'BASIN BOUNDARY POINTS NOT WELL DEFINED') 
RETURN1 

9200 WRITE( IOUNIT, 92 10) 

9210 F0RMAT(/1X, 'NO DATA AVAILABLE') 

RETURN 1 


END 
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C *** SUBROUTINE PCINS *** 

C 

C OBJECTIVES: 

C STOP PROGRAM - NO OTHER INPUT 

C 

C 

SUBROUTINE PCINS( IOUNIT) 


WRITE ( IOUNIT , 1 ) 

1 FORMAT(//15X, '*** END ***') 
C 

RETURN 

END 
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ERRATA 

FOR 

COMBINING REMOTELY SENSED AND OTHER 
MEASUREMENTS FOR HYDROLOGIC AREAL AVERAGES 
NASA-CR170457 OCTOBER 31, 1982 


Page No. 


16 

1st line 

Paramater rather than parameters 

29 

3rd para 
2nd sent 

should read --Advanced very high 
resolution radiometers (AVHRR) . . 

30 

1st para 
last sent 

ADD -- in the visible and 8 km 
in the thermal infrared. 

38 

sect. 5.1 

Change second bullet to read: 
• NOAA-8 AVHRR, and 

39 

2nd para 
1st sent 

Change to read: The NOAA-8AVHRR 
on the TIROS satellites... 

42 

Eq (3 . 1) 

Change last part of right hand 

r T 



to read: ( l-A^-A^) 



60 

Eq (3.5) 

Eliminate second = sign 


64 

2nd para 
2nd line 

put a subscript a with C at end 
of sentence 


73 

3rd para 
3rd sent 

should read: The 17 points along 
the line... 

w 

77 

last para 
2nd sent 

delete the phase -- can be increased 

• 

78 

2nd para 
4th line 
5th, 7th & 
9th lines 

Schmugge not Schmugg 

T's should have subscript B 


81 

2d & 3d para 

The two a ' s should be a ' s 


82 

2d para 
3rd line 

spelling of function 


B-l 


B-l 

1st para 


5th sent 

C-l 

last sent 

D-19 

Eq (D. 33) 

F-17 

Eq ( F . 19 ) 


spelling of approach 
means rather than mean 
should read: 



1st 3 below line should be 
squared (3^) 
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